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Numerical  solutions  for  flow  and  heat  transfer  in  a 
single-pass,  return-flow  heat  exchanger  are  given  for 
Reynolds  numbers  (Re)  of  100,  500,  and  1000,  and  Prandtl 
numbers  (Pr)  of  0.7,  20,  and  100.   The  flow  is  modeled  as 
laminar,  axisymraetric,  incompressible,  and  in  steady  state, 
and  the  Navier-Stok.es  equations  are  expressed  in  terms  of 
stream  function  and  vorticity. 

Prior  to  the  solution  of  the  problem  by  means  of  a 
finite  difference  method,  the  original  geometry  of  the 
system  in  the  physical  plane  is  transformed  into  a 
rectangular  domain  with  square  grids  in  a  computational 
plane.   Grid  clustering  techniques  are  also  employed  to 
reposition  grid  lines  for  good  resolution  in  those  regions 
where  the  flow  and  temperature  fields  change  most 
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rapidly.   The  original  partial  differential  equations  are 
expressed  in  finite-difference  forms,  and  in  order  to 
assure  convergence  in  the  numerical  iteration,  the  matrix 
coefficients  in  the  matrix  equations  are  reconstructed  to 
make  those  diagonal  terms  dominant  in  the  coefficient 
matrix. 

Numerical  data  indicate  that  the  flow  is  accelerated 
in  the  end  cap.   There  is  a  recirculation  zone  near  the 
inner  end  of  the  inner  pipe,  where  the  main  flow  is 
separated  from  the  wall.   The  mean  Nusselt  number  for  the 
exchanger  is  found  to  be  proportional  to  the  Reynolds 
number  and  the  Prandtl  number.   The  present  exchanger  has  a 
mean  Nusselt  number  that  is  about  1.57  to  1.87  times  that 
for  heat  transfer  in  a  singular  pipe  under  similar 
operating  conditions  (Re  and  Pr  number).   The  exchanger  is 
thus  shown  to  be  effective  in  heat  transfer. 
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CHAPTER  I 
INTRODUCTION 


The  single-pass,  return-flow  (SPRF)  heat  exchanger 
shown  in  figure  1  consists  of  a  circular  pipe  inserted 
inside  an  outer  pipe.   The  outer  pipe  is  sealed  at  one  end 
so  that  fluid  entering  the  inner  pipe  may  be  deflected  into 
the  annulus  formed  between  the  inner  and  the  outer  pipe  of 
the  exchanger.   Either  the  inner  or  the  outer  pipe  wall,  or 
both,  may  be  heated,  and  the  heat  input  may  be  simulated  by 
different  conditions  at  these  walls. 

The  heat  exchanger  described  above  is  structurally 
similar  to  the  bayonet  tube  studies  by  Hurd  [1]  in  1946; 
however,  they  serve  different  purposes.   In  the  bayonet 
tube,  the  double-pipe  construction  is  essentially  used  for 
transporting  fluid.   Hence,  the  entering  and  exiting  fluids 
exchange  heat  with  each  other  via  the  heat  conducting  inner 
wall,  and  whether  heat  is  transferred  in  or  out  of  the 
system  depends  on  the  temperature  difference  between  the 
fluid  in  the  annulus  and  the  medium  surrounding  it  at  the 
shell  side.   This  is  not  so,  however,  for  the  heat 
exchanger,  in  which  the  inner  pipe  wall  may  be  heated. 
Then,  the  entering  and  exiting  fluids  both  exchange  hent 
with  this  wall,  and  the  system  becomes  a  compact. 


c 
o 

•H 
•P 
CO 
bO 

•H 

m 

0) 

> 


0) 

T3 

C 

e 
<u 
•p 

CQ 

!>, 
W 

(U 

x: 

W) 

c 

•H 

IS 

o 

Si 
CO 


+> 
CO 
B 
(U 
jC 
O 

CO 


•H 


stand-alone  device.   In  fact,  the  exchanger  possesses  all 
tne  virtues  of  the  bayonet  tube,  such  as  the  freedom  from 
stress-and-strain  problems  in  the  inner  pipe  and  the  ease 
of  replacement  of  tubes  in  case  of  failure,  while  not 
suffering  from  the  lack  of  versatility  that  usually 
handicaps  the  use  of  the  bayonet  tube. 

A  review  of  the  literature  reveals  that  there  is  a 
lack  of  basic  research  dealing  with  the  heat  exchanger. 
Hurd  [1]  studied  the  mean  temperature  difference  in  the 
bayonet  tube  for  four  different  arrangements  of  fluid  flow 
directions.   Large  temperature  differences  were  found  when 
flow  direction  in  the  annulus  was  counter  to  the  flow 
direction  of  the  heat  exchanger  medium  at  the  shell  side. 
A  more  recent  study  dealt  with  the  use  of  the  heat 
exchanger  in  solar  collectors  [2j.   In  this  application, 
the  outside  pipe  of  the  exchanger  becomes  the  inner  wall  of 
a  dewar,  and  the  outside  surface  of  this  pipe  is  coated 
with  a  spectral  selective  coating  to  enhance  its  heat 
absorption.   Owens-Illinois  has  used  this  concept  in 
developing  its  SUNPAK  collectors  with  success  [3]. 

It  is  quite  natural  for  one  who  is  interested  in  fluid 
flow  and  heat  transfer  analysis  to  consider  the  SPRF 
exchanger  as  a  combination  of  two  separate  components, 
circular  pipe  and  annulus,  because  their  fluid  flow  and 
heat  transfer  have  been  thoroughly  studied  in  the 
literature,  e.g.  [4].   However,  because  of  the  presence  of 


the  round  cap  at  one  end  of  the  heat  exchanger,  the  fluid 
flow  in  the  upstream  region  may  be  affected  by  that  in  the 
downstream  region.   As  a  result,  to  divide  the  exchanger 
into  two  subsystems  and  to  patch  up  their  solutions  at 
points  of  matching  may  be  an  oversimplification  of  this 
complicated  problem.   A  realistic  modelling  of  the  system 
will  require  the  treatment  of  the  system  as  a  single 
unit.   This  motivates  the  research  of  this  dissertation. 

In  this  work,  the  flow  is  assumed  to  be  laminar, 
axisyrametric,  incompressible,  and  in  steady  state.   Stream 
function  and  vorticity  are  used  in  the  solution  of  the 
problem  (Chapter  II).   The  original  irregular  geometry  of 
the  system  in  the  physical  plane  is  transformed  to  a 
rectangular  domain  with  square  grids  in  the  computational 
plane.   Grid  l.nes  are  clustered  to  provide  good  resolution 
in  regions  where  flow  and  temperature  fields  change  most 
rapidly  (Chapter  III).   The  problem  is  solved  numerically 
with  an  iterative  scheme  and,  in  order  to  assure  a 
convergent  solution,  diagonal  matrix  elements  in  the 
governing  matrix  equation  are  maximized  (Chapter  IV). 
Pressure  distribution  and  Nusselt  number  are  evaluated  at 
Re=100,  500,  and  1000,  and  Pr=0.7,  20,  and  100  (Chapter 
V).   Thus,  laminar  flow  is  stable  and  the  use  of  the  steady 
flow  form  of  the  governing  equation  is  physically 
credible.   Numerical  results  indicate  that  the  SPRF 
exchanger  is  superior  to  heat  transfer  in  a  circular 


pipe.   The  present  exchanger  has  a  mean  Nusselt  number  that 
is  about  1.57  to  1.87  times  that  for  heat  transfer  in  a 
singular  pipe  under  similar  operating  conditions  (Chapter 
VI). 


CHAPTER  II 
FORMULATION  OF  PROBLEM 


The  system  under  investigation  is  shown  in  figure  2, 
where  the  dimensional  ratios  of  the  pipes  are  also  provided 
in  the  figure.   Note  that  the  thickness  of  the  inner  pipe 
is  only  one-hundredth  of  its  radius,  while  the  inner  end  of 
the  pipe  that  is  facing  the  spherical  cap  is  rounded.   The 
radius  of  the  inner  pipe  was  chosen  such  that  the  cross 
section  of  flow  in  the  pipe  is  equal  to  the  ring  passage 
between  the  pipe  and  the  cap,  and  this  area  is  also  equal 
to  the  cross  section  of  the  annulus  region;  see  dashed  line 
marks  in  the  figure.   The  length  of  the  pipes  is 
approximately  twelve  times  the  radius  of  the  inner  pipe,  a 
length  short  enough  to  show  the  entrance  effect,  yet  long 
enough  so  that  the  boundary  conditions  imposed  at  the  exit 
section  of  the  annulus  will  have  little  effect  on  the  flow 
conditions  inside  the  annulus. 

II . 1   Governing  Equations 
For  the  problem  under  investigation,  the  flow  is 
treated  as  steady,  incompressible,  and  laminar.   The  fluid 
medium  is  assumed  to  have  constant  thermophysical 
properties.   Both  the  flow  and  the  temperature  field  are 
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invariant  in  the  ^   direction;  the  flow  is  thus 
axisymmetrical ,  which  also  permits  the  use  of  the  stream 
function  in  the  analysis.   It  is  further  assumed  that  the 
Eckert  number  is  small  enough  so  that  dissipation  is 
negligible.   Under  these  conditions,  the  momentum  equations 
in  r  and  z  directions  can  be  combined  to  derive  the 
vorticity  transport  equation  [5] 


—  il>  n     -1^0.     +  ^  ^p     =4—  (^  +  ^       +  —  a  -  -^)   (2.1) 
r   r  z   r   z  r    2  z   Re   zz    rr    r   r    2 
r  r 


where  the  vorticity  f2  is  related  to  the  stream  function  ij) 
by  the  following  relationship: 


*^^  +  'i'^^  -  T-  "^t.  =  -  r*"  (2.2) 

zz    rr   r   r 


As  usual,  the  stream  function  is  related  to  the 
velocity  u  and  v  by 


u  =  7  'l'^  (2.3) 

V  =  -  f  1^2  (2.4) 


The  energy  equation  can  be  expressed  in  terms  of 
stream  function  as 


2.  ip      Q      _J_^   0   =  1_  (0    ^.0    +J_0)         (2.5) 
r   r   z   r  ^z   r   Pe  ^  zz    rr    r   r'^         \^'JJ 


where  Pe  stands  for  the  Peclet  number,  a  product  of  the 
Prandtl  number  and  the  Reynolds  number,  Pe=PrxRe. 

In  the  formulation  given  above,  the  spatial  variables 
r  and  z  have  been  made  diraensionless  by  dividing  the 
original  r   and  z   by  rj_,  the  radius  of  the  inner  pipe. 
The  velocities  u  and  v  have  been  normalized  with  reference 
to  the  mean  velocity,  Ujj,;  0  is  a  dimensionless  temperature, 
defined  as  e=(T^-T)/(T^-TQ) ,  where  T^  and  Tg  refer  to  the 
inner  wall  temperature  and  the  fluid  inlet  temperature, 
respectively.   There  are  five  unknowns  {n,v ,Q,^ ,^)    in  these 
equations,  and  there  are  five  equations  to  solve  them. 

I I. 2   Boundary  Conditions 
Because  of  the  large  number  of  conditions  to  be  used 

in  the  analysis,  it  is  convenient  to  summarize  them 

according  to  the  locations  where  they  will  be  applied. 

Inlet:   The  velocity  and  temperature  are  both  uniform  at 
the  inlet.   Here  the  v  component  of  velocity  is 
zero. 

Exit:    Because  of  the  exit  section  being  remote  from  the 
inlet  section,  u,  v,  t|j ,  and  9.   are  considered 
invariant  in  the  z  direction  at  exit.   This 
assumption  has  been  supported  by  a  study  in  which 
it  showed  that  the  exact  nature  of  the  downstream 
condition  had  little  effect  on  the  solution  of  the 
problem  [6].   It  is  further  assumed  that  conduction 
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in  the  fluid  is  small  at  the  exit  section,  a  valid 
assumption  for  Pe  >  100  as  proposed  by  Singh  [?]• 

Axis  of  symmetry:   The  first* derivative  of  u  and  Q  with  r 
must  be  zero  along  the  axis  of  symmetry.   Here  the 
velocity  v  is  also  zero. 

Walls:   Both  u  and  v  velocities  vanish  at  the  inner  and 
outer  pipe  walls.   Physically,  to  satisfy  these 
conditions,  the  vorticity  must  be  continuously 
generated  at  these  walls.   For  heat  transfer 
analysis,  a  constant  temperature  condition  is 
imposed  at  the  inner  wall,  while  the  outer  wall  is 
insulated. 
The  conditions  mentioned  above  can  be  formulated  as 

follows. 

Inlet: 

u  =  1  (2.6a) 

V  =  0  (2.6b) 

0  =  1  (2.6c) 

i>   =  j-  (2.6d) 

n  =  0  (2.6e) 


Exit; 


u   =v  =    ^      =    Q     =0  (2.7a) 
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\z    =0  (2.7b) 


Axis   of   symmetry: 


u      =v    =    e      =    ri,    =    Q.   =   0  (2.8) 


Inner  pipe  wall: 

u=v=0=O  (2.9a) 

4-  =  J  (2.9b) 

Outer  pipe  wall: 

u  =  v  =  G=i(;=0  (2.10) 

n 

where  the  subscript  n  for  0  denotes  the  normal  derivative. 

The  missing  vorticity  conditions  at  the  walls  may  be 
derived  by  using  Eq.  (2.2),  in  which  a  great  simplification 
can  be  made  by  using  the  facts  that  (i)   u  =  v  =  0  at  the 
walls,  and  (ii)   the  walls  may  be  regarded  as  one  of  the 
streamlines.   For  the  time  being,  this  condition  is 
expressed  as 

n=-[  —  (i|;    +   ^        _1^)]  (2.11  ) 

w       r   zz    rr   r  ^r   wall 

It  will  be  simplified  later  in  Chapter  IV. 


CHAPTER  III 
GRID  GENERATION  AND  CLUSTERING 


Because  of  the  irregularity  of  the  geometry  under 
investigation,  the  solution  of  the  problem  in  the  original 
physical  plane  is  difficult.   It  is  desirable  to  use  grid 
generation  techniques  to  transform  the  grid  points  in  the 
physical  plane  to  the  computational  plane  so  that 
computation  can  be  conveniently  carried  out  in  the  new 
plane . 

Refer  to  the  mapping  diagrams  shown  in  figure  3.   The 
physical  plane  is  in  (z,r)  coordinates,  while  the 
computational  plane  is  in  (£;,n)  coordinates.   It  is 
necessary  to  transform  the  domain  abcdefghija  in  the 
physical  plane  to  the  domain  abcdefghija   in  the 
computational  plane,  the  mapping  being  one-to-one.   Notice 
that  the  fluid  enters  the  heat  exchanger  through  the  inner 
pipe,  and  the  fluid  actions  are  more  intense  near  the  wall 
and  around  the  turning  point  (h)  at  the  tip  of  the  inner 
pipe;  a  fine  grid  is  necessary  to  study  the  fluid  flow  and 
heat  transfer  in  these  regions.   To  facilitate  grid 
generation  to  meet  these  requirements,  the  total  system  is 
subdivided  into  three  domains:   domain  A  in  bcdghib,  domain 
B  in  abija,  and  domain  C  in  fgdef;  and  transformations  will 
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Figure  3   A  Schematic  Showing  Mapping  from  the  Physical 
Plane  to  the  Computational  Plane 
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be  made  separately  in  these  regions  to  map  the  grid  points 
from  one  plane  to  another.   Also  note  that,  after 
transformation,  the  streamlines  in  the  physical  plane  will 
correspond  roughly  to  the  constant  n  lines  in  the 
computational  plane. 

The  transformation  will  be  accomplished  in  three 
steps.   In  the  first  step,  the  domain  A  will  be  mapped  to 
domain  A*.   The  next  step  is  to  fine-tune  those  grid  lines 
in  this  domain  so  that  a  dense  grid  is  formed  close  to  the 
walls  and  near  the  turning  point  of  the  inner  pipe. 
Finally,  grids  in  domains  B  and  C  are  filled  in  so  that 
they  match  those  grids  already  laid  out  in  domain  A.   The 
details  of  these  transformations  will  be  spelled  out  in  the 
three  sections  that  follow. 

I II . 1   Grid  Generation  in  Domain  A 
Differential  equation  methods  will  be  used  to  generate 
grids  in  domain  A.   Specifically,  the  elliptic  partial 
differential  equation  technique  developed  by  Thompson, 
Thames,  and  Mastin  (TTM)  will  be  used  [8j.   In  order  to 
simplify  computation  of  the  velocity  and  temperature  fields 
to  be  made  later,  the  grids  in  domain  A   in  the 
computational  plane  are  taken  to  be  squares,  and  Laplace 
equations 


5    +    K        =   0  (3.1) 

zz    rr 
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and 


n    +  n    =0  (3.2) 

'zz    rr  v^"^y 


are  used  in  the  physical  plane  to  generate  those  square 
grids  in  the  computational  plane;  see  figure  3.   It  can  be 
shown  (see  Appendix  A)  that  these  equations  are  transformed 
to  the  following  two  in  the  computational  plane 


"^C    -    2^^.    ^    ^^n    =    °  (3.3) 

ar^^    -   26r^      +    yr        =0  (3.4) 

55  5n  nri 

where 

2  2  2  2- 

a   =   z      +    r      ,    3    =   z   z      +    r    r      ,    y   =   z      +    v  (3.5) 

Tin  i    T]  E.    r\  E,  t, 


It  is  imperative  that  the  physical  boundary  of  the 

■X-   ^   -X- 

inner  pipe,  that  is,  ihg,  be  mapped  onto  i  h  g  ,  and  bed  be 

-K-   -X-   -)(- 

mapped  to  b  c  d  ;  both  are  constant  n  lines  in  the 
computational  plane.   Similar  relations  must  be  established 
for  the  bi  and  dg  lines,  which  become  constant  5  lines  in 
the  computational  plane.   These  one-to-one  mapping 
relations  are  usually  referred  to  in  the  literature  as 
boundary  conditions  for  the  transformation  equations  (3.3) 
and  (3.4).   In  actuality,  they  provide  the  necessary 
conditions  so  that  constant  n  lines  conform  to  the  physical 
boundary  of  the  actual  system. 

To  lay  out  z  coordinates  that  enable  a  dense  grid  to 
be  formed  near  the  turning  point  h,  a  cubic  polynomial 
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equation 


z,  =  Zq  .  a^.,  .  a^tj  .  a^^l  (3.6) 


proposed  by  Nietubicz  et  al.  is  used  [9].   The  Zq  and  Zj_ 
refer  to  the  z  positions  at  index  iQ  and  i^,    respectively; 
Iq   refers  to  the  index  at  the  turning  point.  "V^   is  a 
normalized  index,  given  as 


1 .  -  ip, 
4-.  =  -r^ r^  (3.7) 


1         Ir.        -        1 


0 


Equation  (3.6)  permits  the  z  coordinates  to  be  laid  out 
with  the  size  of  the  increment  Az  increasing  with  the  value 
of  the  index  ij_.   Once  the  positions  of  the  first  and  last 
point  on  the  z  axis  are  fixed,  and  correspondingly,  values 
for  the  indices  Iq  and  i^  are  assigned,  and  with  the 
provision  of  the  sizes  of  the  first  and  the  last  z 
increment  as  additional  input,  this  equation  can  be  used  to 
generate  an  array  of  z  positions  having  spacings  as 
desired.   The  use  of  this  equation  is  discussed  fully  in 
Appendix  B. 

Positions  of  interior  nodal  points  (z,r)  are  found  by 
solving  Eqs.  (3.3)  and  (3.4).   Since  these  equations  are 
nonlinear,  a  numerical  technique  must  be  used  to  solve 
them.   In  this  method,  all  derivatives  in  these  equations 
are  represented  by  central  differences,  and  the  resulting 
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difference  equations  take  the  following  form: 

az._,^.-2(an)^i,j-az,^^^.  =  "Y  (  ^i  ,  j_,  ^^i  ,  j,i  )  ^  f  ^.^  (3-8) 
ari_,,j-2(a^y)r,^..ar,,,^.  =  "T  (  ^i  , -j_,  ^i  ,  j.^  )  -  i  ?^^  (3-9) 


where 


«  =Tt^^i, 0^1-^1.3-1^'^  (^i,d.i-i,d-i)'^     (5.10) 
e  =T  t(^i.i,j-"i-i,j^("i,j^i-"i,j-i^ 


-h^-u^,^^i-^,j^     ^  ^'u^,3-'i-^,3^  ^ 


(3.12) 


^n  ^  ("i^1,j^1-^i-1,j-l'"i-1,j-l"^i-1,j-1^      (^•''^^ 


and 


^5n  =  (^i.1.j.1-^i.1,j-1^^i-1,d-1-^-1,j.l)      (^-'^^ 

Iteration  is  used  in  conjunction  with  successive  line 
over-relaxation  to  solve  these  equations.   Tne  relaxation 
procedure  is  formulated  as  follows  for  each  ti  line  (i.e., 
for  a  fixed  j  index): 


z.  .  =  z.  .  +  w^  (z.  .-z.  .)  (3.15) 


and 


r^"''.  =  r^  ,  ^  w^  (r.  ,-r^  ,)  (3.16) 
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where  1  <  i  <  imcv   In  these  equations,  z-  ■  and  r,-  .:  are 
met  A  -'-  >  J        ■'-  f  o 

the  solutions  of  Eqs.  (3.8)  and  (3.9);  w^  and  Wj,  are 
relaxation  parameters  for  z  find  r.   The  superscript  k  for  z 
and  r  refers  to  the  Iteration  count.   In  the  iteration,  the 
initial  guess  for  the  grid  locations  is  made  by  using 
linear  interpolation  of  those  grid  already  located  at  the 
boundaries.   It  is  found  that  a  relaxation  parameter  of  1.7 
yields  fast  convergence  for  both  w^.  and  Wj, ,  and  the 
iteration  is  terminated  if 


k+1   k 
z .  .-z . 


<  10 


-4 


and 


k+1   k 
r .  .-r .  . 
1,3   i.J 


<  10 


-4 


The  grid  generated  by  the  TTM  solver  is  shown  in 
figure  4  (a).   The  spacing  along  the  z  axis  is 
satisfactory;  however,  the  spacing  in  the  r  direction  is 
not,  which  will  be  corrected  later.   It  should  be  repeated 
that  these  grids  are  generated  by  using  square  grids  in  the 
computational  plane.   Also  notice  that,  for  the 
transformation  made  that  is  illustrated  in  this  figure,  103 
points  are  used  along  the  z  direction,  41  points  in  the  r 
direction. 


III. 2   Grid  Clustering  in  Domain  A 
The  grids  generated  in  the  preceding  section  use 
Laplace  equations  for  transformation.   If  Poisson's 
equations  are  used  instead,  the  source  terms  in  the 
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Poisson's  equations  may  be  used  to  control  grid  positions, 
then  a  separate  grid-line  clustering  effort  may  be 
spared.   Sorenson  and  Steger* found  out  that  this  method  of 
clustering  grid  lines  may  lead  to  complexities  or  even 
unsatisfactory  results  in  the  computation  [10].   For  this 
reason,  the  grids  are  not  transformed  using  Poisson's 
equations  in  this  work;  Laplace  equations  are  used  instead, 
and  the  generated  grid  lines  are  clustered  later  by  means 
of  clustering  functions.   In  this  effor",  Sorenson  and 
Steger's  method  will  be  used  as  described  below. 

The  grid  lines  in  the  z   direction  have  been  laid  out 
satisfactorily;  no  changes  need  be  made  of  them.   Lines  in 
the  r  direction  must  be  clustered  so  that  a  dense  grid  is 
where  the  change  in  the  gradient  of  the  flow  variable  is 
most  rapid.   In  other  words,  lines  of  constant  n  in  the 
physical  plane  must  be  repositioned  for  each  line  of 
constant  C   This  clustering  of  n  lines  can  be  accomplished 
in  two  steps  as  follows. 

(1)   The  arc  length  along  each  constant  5  line  in 
figure  4  (a)  is  computed  with  the  following  relation: 


S.  .  ,  =  S.  .+[(z.  .  ,-z.  .)^+(r.  .  ,-r.  .)^]^^^  (3.17) 

i,J+1     i,J     1,0+1   1,0      1,0+1   1,0 


Since  the  first  point  is  indexed  1  ,  Sj_  -,  =0  for  the  sake 
of  consistency.   Note  that  j  refers  to  the  index  for  the 
jth  line  of  constant  n ;  the  arc  length  is  computed  for  a 
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fixed  i  (or  constant  5).   Equation  (3.17)  can  be  used  to 
compute  the  arc  length  for  each  grid  point  in  figure  4  (a), 
and  this  yields  a  table  of  z,-  .;  and  r,-  a   as  functions  S-;  ^. 

This  part  of  the  work  is  essentially  used  for  locating 
grid  positions  in  the  physical  plane.   Once  they  are  found, 
the  constant  ri  lines  in  the  physical  plane  may  be 
discarded. 

(2)   The  new  grid  lines  in  the  physical  plane  will  be 
laid  out.   An  exponential  stretching  function  is  used  so 
that  the  grid  spacing  between  constant  n  lines  increases 
monotonically  with  the  value  of  j.   The  following  relation 
meets  this  requirement. 


h.i.^   =  s^_  .  .  aSqCUcjJ-''  (3.18) 


Once  again  S.  .  =  0.   In  this  equation,  aS„  is  the  desired 

minimum  grid  spacing  near  the  wall  of  the  outer  pipe;  e  is 

chosen  such  that  S.  .     =  S.  .    .   Equation  (3.18) 

l"i  in  '.  -'         I 

'■^max  ''-'max 

enables  arc  lengths,  along  a  constant  5  line,  to  be  found 

that  give  the  desired  grid  locations  for  new  ^  lines.   The 

actual  positions  for  these  grids  in  the  physical  plane  are 

found  oy  interpolating  the  Zj_  ^  and  rj_  ^  tables  constructed 

earlier . 

There  is  one  parameter  £  that  remains  to  be  found  in 

Eq.  (3.18).   As  has  been  mentioned,  this  parameter  controls 

the  position  of  the  outermost  boundary,  so 
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F(e)  =  S.  .     -  S        =0  (3.19) 

■^"^raax      '"'max 


The  Newton-Raphson  method  may  be  used  to  find  this  e.   In 
this  method,  the  nth  iterated  e  is  expressed  as 


F(e"~^) 
,n  ^  ^n-1 3^_^  (3.20) 

^    ^     F'(sr') 


1 


where 


FC."--')  =  S,  ,     .  ^  UU^I-h^"'^'' -U  (3.21 


and 

(e.   ) 


X[ernj._-2)-1]}  (3.22) 


1    "max 


Notice  that  this  e^   must  be  found  for  each  E.    line  in  the 
reconstruction  of  new  grids. 

This  clustering  scheme  described  above  monotonically 
increases  the  grid  spacing  away  from  the  outer  pipe  wall. 
This  creates  a  slight  problem  for  those  grids  near  the 
inner  pipe  wall  because  the  rate  of  change  of  the  gradient 
of  the  flow  variable  is  also  rapid  there.   It  is  necessary 
to  cluster  the  grid  lines  near  boundary  ihg  (figure  3)  for 
good  resolution.   To  accomplish  this,  an  image  method  is 
employed  as  follows. 
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In  the  image  method,  the  S.  .     in  Eq.  (3.18)  is 

'^max 
redefined  to  be  S.  .    /2,  which  length  corresponds  to 

'•^max 
that  length  to  index  ( Jniax"^^  ^Z^'  ^"^  ^°    that  line  at  its 

center.   Then  lines  [  ( Jmax'^'' ^/^"^^  ^  ^°  Jmax  "^^d  not  be 

regenerated  but  laid  out  by  regarding  them  to  be  the  mirror 

image  of  those  lines  from  [{j^g^^  +  '])/2-']]    to  1. 

The  clustered  grid  lines  are  shown  in  figure  4  (b), 

where  ASq  is  0.25^  of  the  inner  pipe  radius.   An  enlarged 

grid  near  the  turning  point  is  shown  in  figure  5.   The 

clustered  grid  lines  appear  to  be  uniformly  distributed 

near  the  pipe  wall  and  the  turning  point  of  the  pipe.   The 

grid  generation  task  in  domain  A  is  thus  completed. 

III. 3   Grid  Generation  in  Domain  B  and  C 
The  grid  generation  in  domain  B  and  C  is  rather 
straight-forward  once  the  grids  are  laid  out  in  domain  A. 
Recall  that  the  positions  of  the  grid  lines  in  the  z 
direction  in  the  physical  plane  are  generated  by  using  Eq. 
(3.6),  in  which  values  for  Zq,  z^,  Iq,  i^,  Az^ ,  and  Az^ 
must  be  provided  as  input;  the  grid  spacing  increases  with 
the  index  i.   Then,  in  domain  B,  the  z  positions  for  the 
grid  lines  must  be  laid  out  from  a  to  b  because  close  to 
the  inlet  section  a  dense  grid  is  desired.   In  domain  C, 
however,  the  direction  of  z  must  be  reversed  because,  at 
the  exit  section,  the  flow  should  approach  fully  developed; 
a  dense  grid  there  is  unnecessary. 
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As  for  the  lines  of  constant  n  in  domains  B  and  C,  a 
straight  ray  extension  of  those  in  domain  A  is  adequate  for 
this  study.   This  method  of  grid  generation  results  in  an 
orthogonal  non-uniform  grid  net  in  domains  B  and  C.   A  plot 
of  the  finished  grid  is  shown  in  figure  6.   Because  domains 
B  and  C  are  now  linked  to  domain  A,  there  are  a  total  of 
163  points  in  the  z  direction;  the  number  of  points  across 
the  flow  passage  in  the  r  direction  is  still  41,  which  is 
unchanged. 

III. 4   Evaluation  of  the  Source  Terms 
in  Poisson's  Equations 

The  dense  grids  near  the  walls  and  the  center  line  of 
the  exchanger  may  be  considered  a  result  of  heat  sources 
and  sinks.   Hence,  although  Poisson's  equations  were  not 
used  to  generate  grids,  the  redistributed  grid  locations 
implicate  that  these  source  terms  do  present  now;  the 
original  Laplace  transformation  relations  are  no  longer 
valid,  and  the  forcing  functions  in  the  Poisson's  equations 
must  now  be  determined. 

It  has  been  shown  in  Appendix  A  that  these  Poisson's 
equations  in  the  computational  plane  are  of  the  following 
form: 


az_.  -  26z.   +  TZ^^  =  -  J^(Pz^+Qz  )         (3.23) 
55      Cn     nn  S    n 


ar,,  -  2Br,   +  yr    =  -  J^(Pr,+Qr  )         (3.24) 
55      5n     nn  <;    n 
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where  J  is  the  Jacobian  of  transformation 


^   =   ^5^  -  ^n^? 


(3.25) 


Equations  (3.2:))  and  (3.24)  may  be  used  to  derive  the 
source  expressions  as 


v/here 


and 


P  =  (z^Dr  -  r^Dz)/J- 


=  (rpDr  -  z^Dz)/J- 


Dr  =  ar,^  -  26r^   +  yr 

CC       Cn      r\r] 


Dz  =  az^^  -  26z^   +  YZ 

5?      5n      nn 


(3.26) 
(3.27) 


Equations  (3.26)  and  (3.27)  are  to  be  solved  numerically  to 
find  the  forcing  terras. 


CHAPTER  IV 
NUMERICAL  SOLUTION 


With  the  grids  generated  in  the  physical  plane  (z,r), 
it  is  only  necessary  to  solve  the  problem  in  the 
computational  plane  (5,n)  so  that  a  transformation  back  to 
the  physical  plane  gives  the  final  answers  to  the 
problem.   The  numerical  solution  of  the  problem  is  provided 
in  this  chapter.   The  governing  equations  and  boundary 
conditions  in  the  transformed  plane  will  be  given  first  in 
the  section  that  follows. 


IV. 1   Transformed  Governing  Equations 
and  Boundary  Conditions 

The  governing  equations  in  the  computational  plane  are 

derived  by  using  the  partial  differential  identities  given 

in  Appendix  A. 

Vorticity    transport   equation: 

t2  5        Jz  p        Jz. 

-^^W-^,\^    .^(.^,^-.^,^)]  (4.1) 
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Vorticity  definition  equation: 

Jz 


Jz 


c^^^^    -    23^^^  .  Y^^^  .(J^P.  ^)*^  .(jV  -^)^^ 


=    -    J    VQ. 


(4.2) 


Velocity  equations 


1 


^  =  jF  ^^5^  -  ^^^ 


(4.3) 
(4.4) 


Energy  equation: 

«^?  -  230^^  .  ye^^  .(j2p-  4^)0^  .(j2q.  _!l)e^ 


i^(,^e^  -^%) 


(4.5) 


Boundary  conditions  can  be  derived  accordingly 
(Appendix  A).   Notice  that,  in  the  application  of  boundary 
conditions,  all  conditions  of  the  first  kind  will  be  used 
as  is;  no  transformations  of  them  need  be  made.   Other 
conditions  must  be  changed  to  the  computational  plane,  and 
they  are  listed  together  with  those  first-kind  conditions 
below. 


Inlet; 


u  =  1 
V  =  0 
0  =  1 


(4.6a) 
(4.6b) 
(4.6c) 


^  =   2 
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li  (4.6d) 


Exit; 


a  =  0  M  (4.6e) 


u   =v  =   i>      =   a      =  0  (4.7a) 

9   =  -^  e  (4.7b) 


Axis  of  symmetry; 


z  u  -  z  u,  =  0  (4.8a) 

V  =  0  (4.8b) 

z  0   -  z  0   =0  (4.8c) 

^],  =   a  =  0  (4.8d) 


(4.9a) 
(4.9b) 

(4.10a) 
(4.10b) 

(4.10c) 


Inner  pipe 

wall: 

u 

=  V  =  0 

=  0 

'I' 

1 

"  2 

Outer  pipe 

wall: 

u 

=  V  =  0 

Y0   -  60^ 

=  0 

'J' 

-  0 

Notice  that  at  both  inner  and  outer  walls,  the 
vorticity  boundary  conditions  can  be  expressed  as 
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"w  =  --3-^n  ^4.11) 

J  r 


Derivations  of  these  expressions  are  given  in   le  appendix. 

IV. 2   Discretization  of  Governing  Equations 
and  Boundary  Conditions 

Those  governing  equations  given  above  will  be  solved 

numerically.   In  the  solution,  the  grids  are  indexed  as 

(i,j),  with  i  pointed  in  C  direction  and  j  in  n 

direction.   The  governing  equations  are  discretized  as 

follows. 

Vorticity  transport  equation: 

[cy+ReCciy'F.-cisiF^)]^.  .  = 

(C8-ReCl6^^)^^^^^j  +  (C9+ReCl6l'^)^._^^j 
+  {C^0+ReC^6i>.)^.     .    .+(C11-ReCl6ij;.)f^.  .  ^+COn     (4.12) 


Vorticity  definition  equation: 


C1t,  .  .   C2*l.1,o*"ti.i_j.C4*,^j,,.C5*._j., 


+C6fi.  .+CO^r^  (4.1J)) 


Velocity   equations 


u.     .    =   C12^„    -   C13K  (4.14) 

V.     .    =   C14^^    -   C15^g  (4.15) 
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Energy  equation: 

Cie.  .  =  (C8-PeCl6?^)e     +(C9+PeCi6;F_)e 

+  (C10H-PeCl6;F.)e.  .  .  +  (C11-PeCl6^.)e 

+C0  S^^  (4.16) 

In  these  equations,  CO  through  C18  represent  metric 
coefficients  listed  in  Appendix  C.   Those  'I',  ^,    and  Q 
topped  with  bars  are  compact  expressions  for  their 
increments  in  i  and  j  directions;  for  example, 


h  =  h.^,j  -  n-i,j 


and 


*n    "^1,0  +  1    "^1,0-1 


'''^n  -  '''i  +  1,j  +  1  "  *i  +  1,j-1  "■  *i-1,j-1    '^i-1,j  +  1 


It  is  strictly  a  matter  of  computer  programming  expediency 
that  these  substitutions  are  made. 

The  discretized  equations  given  above  can  be  solved  by 
an  iterative  method,  in  which  the  problems  encountered 
include  the  nonlinearity  of  the  vorticity  transport 
equation  and  the  high  convection  rate  in  the  energy 
equation.   Although  use  of  an  under-ralaxation  scheme  in 
the  iterative  method  may  ease  these  problems,  a  small 
relaxation  factor  means  an  excessive  computer  time.   For 
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these  reasons,  the  matrix  coefficients  in  the  matrix 
equations  are  reconstructed  to  make  those  diagonal  elements 
dominant  in  the  coefficient  matrix  so  that  an  iteration 
will  yield  a  convergent  solution. 

The  method  used  for  maximizing  the  matrix  elements  is 
somewhat  similar  to  the  scheme  developed  by  Takemitsu  [11] 
and  is  rooted  in  the  up-wind  scheme  that  is  commonly  used 
in  finite-differencing  flow  terras  in  the  solution  of  fluid 
dynamics  problems  [12].   The  vorticity  transport  equation 
is  expressed  as 


C19^.  .  =  C20f2.  .  .+C21f^.  .  .+C22J^.  .  . 
i,J       1+1 , J     1-1 ,J      i,J+1 

+C23^^^j_^+C0n^^        ■  (4.17) 


This  equation  is  a  compact  version  of  the  vorticity 
equation  derived  in  Appendix  C,  Eq.  (C.35). 
The  energy  equation  is  derived  as 


C240.  .  =  C25e.  .  .+C26e.  .  .+C27Q.  ■  . 
1. J       1+1 ,3  1-1 ,3  i,J+1 

+C280,   .  .^COO^  (4.18) 


In  both  equations,  C19  through  C28  represent 


C19  =  C7+Re[C17K  -  CIS'J^^  +2(|C16'|'^|  +  |ci6;F^|)] 


C20  =  C8-Re(Cl6iF^  -   Cl6ii^^|  ) 


C21  =  C9+Re(Cl6^^  +  Cl64i J  ) 


C22  =  C10+Re(Cl6^^  +   Cl64iJ  ) 
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C23  =  C11-Re(Cl6^^  - 


C16^, 


and 


C24  =  C1+2Pe(Cl6^^   +  |C16^^  ) 


C25  =  C8-Pe(Cl6?^  -  |C16^^  ) 


C26  =  C9+Pe(Cl6<'^  +   C16^^|  ) 


C27  =  C10+Pe(Cl64'.  +   C16^^  ) 


C28  =  C11-Pe(Cl6ii)^  -   C164;^|  ) 


Successive-over-relaxation  (SOR)  methods  will  be  used 
to  solve  these  equations.   In  this  method,  the  SOR  versions 
of  (4.13),  (4.17),  and  (4.18)  are,  respectively, 


'i:l  =  ^^-V^i,o^cf^^<i,o^^^^i",a 


^C4^i,d.1  ^  ^5^i:3-1  ^  ^^'T,l    ^   Kn) 


(4.19) 


1.0 


(^-.)"i,d  ^cw  ^^20.^.1, j  ^  ^21.^:]^. 


+  022^2^^  .  .  +  C23^^^''.  .    +   C0<|^*  ) 
1,0+1        i,J-1       ?^ 


(4.20) 


and 
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*=279j,j,1  *  C28ej;].^  .  COe;^)  (4.21) 


where  w^,  v^,    and  Wq  denote  the  relaxation  factors  for  4), 

f^,  and  0,  respectively.   Superscript  k  represents  the  inner 

iteration  count  in  the  numerical  solution;  the  superscript 

n  for  J^,-  A    in  (4.19)  stands  for  the  outer 
-■■  >  J 

iteration  count  discussed  later.   The  cross-derivative 

_*    _*       _* 
terms  ^ c    >    ^r    t    ai^d  ©^   are  evaluated  using 


T*   _  ,k  k+1        k+1     _  k 

'''en  -  '^i  +  l.j  +  l  "  '^i+Lj-l  ■"  "^i-l.j-l    *i-1,j  +  1 


^^     =   ^  •    .^    •   ^    -   ^  •    .^    •   -,    +   ^  •    -1    ■        y,         -        ^  ■         y,  ■        y, 

Cn    1+1, j  +  1     1  +  1, j_i     i_i,j_i     1-1, j  +  1 


and 


■"^n  -  ^i  +  1,j  +  1  '  ^i  +  1,j-1  "■  ^i-1,j-1  "  ^i-1,j  +  1 


Iterations  sweeps  run  from  i  =  1  to  ijjj^^  and  j  =  1  to  jjjjax* 

Discretization  of  the  governing  equations  is  now 
complete;  attention  will  now  be  directed  to  the  boundary 
conditions  which  must  be  processed  accordingly.   In  the 
discretization  of  boundary  conditions,  it  must  be  noted 
that,  in  the  numerical  solution  to  be  given  later,  ^   and  ^ 
must  be  solved  prior  to  the  solution  of  the  velocity  and 
temperature  fields.   It  is  thus  appropriate  to  discretize 
the  ^   and  ^   conditions  first.   In  doing  so,  conditions  of 
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the  first  kind  will  again  be  used  as  is.   Therefore,  only 
conditions  of  the  second  kind  will  be  addressed.   These 
second-kind  conditions  will  be  discretized  in  the  order  of 
their  appearance  in  conditions  (4.6)  to  (4.11)  below. 

The  exit  conditions  Eq.  (4.7a)  are  discretized  by 
using  a  second-order  central-difference  scheme.   At  the 
exit  section,  for  example,  i|)^  =  0  is  discretized  to  be 

il;.     .     .    =   ^-  ,     ■  (4.22) 

^^max^^'J    '^max-^'J 

This  identity  will  be  used  to  substitute  the  values  of  \li   at 
index  (imax+'''j^  ^^^^   ^^^   difference  equation  (4.13)  is 
evaluated  at  the  exit  section.   This  method  also  applies  to 
n,  u,  and  V. 

The  wall  vorticity  conditions  (4.11)  are  approximated 
by  the  relation 


Here  the  bracketed  terms  are  used  to  approximate  -'|>„^>  s^i^d 
subscript  (w+1 )  for  ^    refers  to  the  grid  location  being  one 
point  away  from  the  wall.   This  equation  has  only 
first-order  accuracy  but  is  adequate  for  use  in  the  present 
work  because  of  its  conservation  property  as  reported  by 
Parraentier  and  Torrance  [13]. 
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Conditions  for  e  and  u  are  treated  next,  and  it  is 
noted  that,  at  the  exit,  the  grids  are  orthogonal  and 
ij;_  _  Q.   Equation  (4.5)  reduces  to 


Jz 

"^5  '   '\n    '  '^''\    '    ^'^^   '   -r-^^  =  ^   %^     ^4.24) 


With  the  use  of  Eq.  (4.7b)  and  after  simplification,  Eq. 
(4.24)  becomes 


Jz 
ye   +  (J^Q  +  — r^)e  =^^  e,  (4.25) 


Central  difference  formula  is  then  used  to  approximate  all 
n-derivative  quantities,  and  0^.  is  discretized  by  the 
second-order  backward  difference,  giving 


G   =^  (Oi    _2  .  -  40.    _^  .  +  30.     .)         (4.26) 
^    ^    ^max  "^'^     ^max  '  '^     ^max'^ 


Along  the  axis  of  symmetry,  Eq.  (4.8a)is  approximated 

by 


^.1  =    (zyz^).1  tUi,2  ^  ^%/^)^-1.1^  ^4.27) 


in  which  the  original  u^  and  u   have  been  approximated  by 


^^)i,1  =  ^.1  -  ^-1,1 
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(u  ).  .  =  u.    -  u.  . 
n  1 ,  I    1 ,^    1,1 

again  a  first-order  approximation. 

Temperature  condition  along  the  same  axis  is  treated 
in  a  similar  manner, 


Qi.1  =  (z  /z  ).1  tei,2  ^  (%/^)Qi-1,lJ  (4.28) 


At   the    outer   pipe   wall, 


^i,1    =Ti7TnT  ^®i,2   ^    (e/^)9i-1,1^  (4.29) 


IV. 3   Numerical  Solution  Procedure 


A  close  examination  of  Eqs.  (4.13)  to  (4.15),  (4.17), 
(4.18),  and  (4.23)  reveals  that 

(1)  the  ^  equation,  (4.13),  contains  all  ^i  but  one  q 
terra; 

(2)  the  Q   equation,  (4. 17),  contains  both  4,  and  q 
terms; 

(3)  ill   and  ^   are  also  related  in  the  wall  vorticity 
condition,  Eq.  (4.23); 

(4)  u  and  v  are  functions  of  ijj  only;  see  Eqs.  (4.14) 
and  (4.15); 

(5)  the  0  equation  (4. 18),  contains  both  4;  and  0 
terras . 
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These  characteristics  implicate  that 

(1)  4)  and  Q.    should  be  solved  for  simultaneously  first; 

(2)  6,  u,  and  v  are  solved  next,  with  the  use  of  ^    as 
input. 

Based  on  these  observations,  an  algorithm  is  developed 
as  follows: 

(1)  make  an  initial  guess  of  i|^  and  Q   values  at  all 
grid  points  and  assign  these  values  as  ^^    a   and 

"Id- 

(2)  update  n  values  at  the  wall  using  (4.23) 


n      2y   /  ,n 
^,w  =-^  ^"^^ 


,  n   s 
j2^  ^^i,w+1  "  '^i,w'' 


(4.30) 


(3)   calculate  new  values  of  ^    from  (4.20)  using  SOR 

until  convergence  is  attained;  assign  these  values 
as  n^  .  and  then  damp  them  by  using 


i»0     1,0 


w, (n .  .  -  n .    . ) 


(4.31) 


(4)   calculate  new  values  of  ij;  from  (4.19),  again  usin^ 
SOR,  until  convergence  is  attained;  assign  these 
values  as  ^^    ^   and  then  relax  them  by  using 


ip.  .    =   ^ .     .+w^(i^.  .    -   ^ .     .) 


(4.32) 


(5)   repeat  steps  (2)  to  (4)  until  both  n   and  ^ 
converge . 
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Notice  that  in  steps  (3)  and  (4)  above,  iteration  is 
used  to  find  convergent  a  and  ^   values  for  each  set  of 
assumed  (or  revised)  and  updated  n  and  ^    values;  these  two 
steps  are  then  called  the  inner  iteration  in  the  numerical 
solution.   Steps  (2)  and  (5),  which  repeat  the  total  effort 
by  renewing  n  and  ^^)   values,  are  thus  called  the  outer 
iteration. 

Based  on  the  experience  of  numerical  work  accumulated 
in  the  solution  of  this  problem,  w,  =  1.5  and  w^  =  0.17 
yield  convergent  solutions;  w^  and  Wo  are  both  taken  to  be 
0.3,  following  the  work  by  Rigal  [14].   These  relaxation 
parameters  are  used  consistently  for  all  Reynolds  numbers 
tested  in  this  work.   In  the  numerical  solution,  results 
for  ij;  and  a   calculated  for  low  Reynolds-number  cases  are 
used  as  the  initial  guess  for  high  Reynolds-number  cases. 
Thus,  for  example,  the  output  of  ^   and  fi  for  Re  =  100  is  used 
as  input  for  Re=500  case  in  step  (1).   For  Re=100,  the 
initial  guess  of  the  vorticity  values  is  taken  to  be  zero 
throughout.   As  for  the  stream  function,  a  uniform  flow  is 
assumed. 

The  convergence  criteria  adopted  are  delineated  as 
follows. 

(i)   The  convergence  for  outer  iteration  is  set  at 


n  +  1    n 
^i,w  -  "i,w 


<  10 
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which  also  enables  the  following  conditions  to  be 
met  [15]: 


,n+1    ,n 

1,J      1»J 


<  10 


and 


1,0     i,J 


<  10 


-4 


(ii)  The  convergence  for  inner  iteration  of  Eqs.  (4.19) 
and  (4.20)  is  set  at 


,k+1    k 

1,0      1,0 


<  10 


-3 


and 


1,0     1,0 


<  10 


-3 


The  inner  iteration  steps  are  limited  to  50  in  order  to 
prevent  it  from  converging  to  an  erroneous  solution  [5]. 
It  was  found  that  once  the  outer  iteration  has  converged  as 
desired,  the  inner  convergence  criterion  can  be  met  with 
just  one  iteration.   A  flow  chart  for  the  numerical 
procedure  is  given  in  figure  7. 

Once  the  stream  function  is  found,  the  velocity  field 
is  calculated  from  Eqs.  (4.14)  and  (4.15)  and  the  energy 
equation  can  be  solved  in  a  straight-forward  manner.   In 
the  solution  of  the  energy  equation,  the  second-kind 
boundary  conditions  are  treated  explicitly.   Miyakoda's 
method  [16]  is  used  in  this  effort,  in  which  the  derivative 
boundary  conditions  are  incorporated  into  nodal  equations 
for  interior  nodes  that  are  located  one  point  away  from  the 
insulated  boundary  (that  is,  the  center  line  or  the  outer 
wall).   Hence,  one  point  away  from  the  center  line, 


42 


Set  initial   guess  fpr  ipl  j  and  ^.  j 


Update  wall   vorticity  Q-  ^^ 


t.; 


Solve  for  vorticity  fi.    .•  by  SOR 


u-n 


yes(^.^j) 


Damped  by  .i;^U^J^j+w^(^i,j-"?,j) 


jk 


Solve  for  stream  function  i|j^.  ^^  by  SOR 


-no- 


yes(i|..^j) 


Damped  by  ^^^5=^;^j+W2(^,.^j-4^';,j) 


I 


■  no- 


Figure  7  Block  Diagram  Showing  the  Numerical  Procedure 
Used  to  Solve  Stream  Function  and  Vorticity 
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tC23  -  (z^/z^).l^Qi,2  =  C24e,^^^2^C25G,_^^2^C260,^3 

---.  -    •-^^7(/;,\,  ^-1,1-CQ^.    (4.33) 


and  one  point  awiay.  from  the-  outer-  wall 


t^23--X67^]ei,2  =  C24e.^^  ^2^0250.  _^^2^C26e.  ^3 

^^27j^f^0i_,,,.COe^^      (4.34) 


It  is  also  noted  that,  before  applying  these  equations 
to  find  0j_  2>  ®i-1  1  must. be  updated  by  using  (4.28)  and 
(4.29).   The  relaxation  parameter  Wq  is  taken  to  be  0.7  for 
all  cases  studied.   A  fast  convergence  results  if  the 
iteration  proceeds  from  i  =  1  to  imax'-  ^^^^   is  following  the 
flow  direction,  and  0  =  Jiiiax  ^°   '' '  that  is,  from  the  inner 
wall  to  the  outer  wall. • 


CHAPTER  V 
EVALUATION  OF  PRESSURE  DISTRIBUTION  AND  NUSSELT  NUMBER 


Of  the  numerous  quantities  of  interest  in  fluid  flow 
and  heat  transfer,  the  pressure  drop  and  the  heat  transfer 
coefficient  are  most  important.   The  flow  and  temperature 
fields  evaluated  in  the  numerical  solution  may  be  used  to 
compute  them  as  will  be  shown  in  this  chapter. 

V.1   Pressure  Distribution 


Pressure  drop  in  the  heat  exchanger  can  be  computed 
with  the  use  of  the  velocity  field  obtained  in  the 
numerical  solution.   However,  because  of  the  geometric 
irregularity  of  the  exchanger,  there  is  a  lack  of  a  common 
axis  along  which  this  pressure  drop  can  be  e.aluated  all 
the  way  from  in^et  to  exit.   xhe  centerline  of  the  pipes, 
which  is  commonly  used  as  a  basis  for  evaluating  the 
pressure  drop  in  pipe-flow  problems,  is  not  suited  for  use 
in  the  present  investigation  because  it  terminates  at  the 
center  of  the  end  cap.   The  inner  pipe  wall,  fortunately, 
has  a  continuous  contour  all  the  way  from  the  inlet  to 
exit;  it  is  thus  used  as  a  basis  for  evaluating  the 
pressure  drop. 
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The  Navier-Stokes  equations  can  be  used  together  with 
the  continuity  equation  and  vorticity  to  derive  the  diraen- 
sionless  pressure  gradients  at  the  wall  as  (Appendix  D) 


8z     Re^r   8r^  ^^'  '  ^ 


and 


where 


8r  ~  Re  dz  ^-^   ■ 


p  =  p*/(pu^^)  (5.3) 


Since  only  the  total  pressure  drop  and  the  pressure 
gradient  along  the  wall  are  of  concern,  pressure  at  exit  is 
assigned  to  be  zero, 


P,^  =  0  (5.4) 


This  eliminates  the  need  for  determination  of  an  arbitrary 
constant  when  the  pressure  drop  is  found  later. 

A  refinement  is  made  by  defining  a  new  diraensionless 
pressure  so  that  the  flow  characteristics  of  both 
Poiseuille  flow  and  annular  flow  may  be  incorporated. 
Notice  that,  for  a  fully  developed  flow  in  an  annulus  [17] 


«^   '""  =  -1  (5.5) 


(8/C*)   dz 
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In  this  equation, 


C^  =  r2[(i+r^2)  _.Jlz£l_]  (5.5) 


°  ln(1/r^) 


and  p^,  is  diraensionless,  p^,  =  p^,/(pu   );  Re  is  the 
Reynolds  number;  r   in  Eq.  (5.6)  is  a  radius  ratio, 

^   =  '^io/'^o- 

For  a  Poiseuille  flow  under  a  fully  developed 

condition,  Eq.  (5.5)  is  changed  to 

Re  ^Pfd  _       ■  f         . 

which  can  be  recast  as 


""^    *"  =  -C*  (5.8) 


(8/C"^)   dz 

Comparing  Eqs.  (5.5)  and  (5.8)  provides  a  clue  for  a 
new  pressure  defined  as 


"^   P  (5.9) 


(S/C"") 

With  this  definition,  the  pressure  gradients  at  the  inner 
wall  become 


and 


12. 
8Z 


a/c"" 


r   8  r 


47 


(5.10) 


12. 


9a 


ar     (8/C'')  8z 


(5.11) 


Correspondingly,  the  exit  pressure  condition  should  satisfy 


p"   =  0 


(5.12) 


following  Eq.  (5.4). 

It  is  also  noted  that,  in  the  solution  of  the  flow 
field,  the  flow  quantities  have  been  assumed  to  be 
invariant  in  the  z  direction.   It  follows  that 


dPex    ^Pfd 


dz     dz 


=  -1 


(5.13) 


must  be  met.  :. 

The  formulation  given  above  completes  the  preparation 
of  evaluating  the  pressure  drop  whose  derivation  will  now 
follow. 

Along  the  inner  pipe  wall  of  the  exchanger,  the 
pressure  drop  can  be  evaluated  using  a  line  integral 


t 


+       +       +      f  <i   dp    J4. 

AP   =  P2  -  P^  =  /t^  -^   dt 


(5.14) 
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where  t  refers  to  a  surface  coordinate.   The  integrand  in 
this  equation  can  be  expressed  by  means  of  a  vector 

0 

analysis  as 


if=vp^  .  t=  (||1;,  .-If  a^)  •  t  (5.15) 


where  t  is  a  unit  vector  tangential  to  the  wall.   Expres- 
sions for  ap'^/az  and  3p"^/9r  have  been  given  previously  as 
Eqs.  (5.10)  and  (5.11).   In  the  transformed  plane,  they 
become  (Appendix  A) 

+       .     .  zivn)    -z  (rf2)^ 
12_  ^  ___[__  1  _l D D L  (5.16) 

3z      (8/C^)  r        J 
and 

+     .     r  Q,  -  r^n 
IS-  =  1 n-S L_IL  (5.17) 

3r    (8/C'')      J 

Notice  that,  in  the  transformed  plane,  the  inner  wall 
corresponds  to  lines  of  constant  n •   The  unit  tangential 
vector  t  in  Eq.  (5.15)  thus  takes  the  following  form 
(Appendix  A): 

A  A 

z^e^  +  r_e^ 
t  =   ^  ^    ^   ^  (5.18) 

;[  y 

The  differential  arc  length  along  lines  of  constant  n  is 
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dt  =   ^  Y   d  ?  (5.19) 

Introducing  Eqs.  (5.15)  through  (5.19)  into  (5-14) 


gives 


.  +     +     + 
Ap   =  P2  -  P-]  = 


^1       (8/C^)         "^"^         ' 


^      "  ^   ^  "-Jr  1  dS  (5.20) 


(8/C^)      J 
which  is  to  be  integrated  by  using  the  trapezoidal  rule. 

V.2   Nusselt  Number 


The  local  Nusselt  number  along  the  inner  wall  can  be 
evaluated  using 


m  m 


where  Cy^   is  a  dimensionless  hydraulic  diameter,  defined  as 


C,  =  D,  /r.  (5.22) 

h    h   1 


Dj^  is  the  hydraulic  diameter;  '^^   in  Eq.  (5.21)  denotes  a 
.nixing  cup  temperature,  defined  as  [18] 
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/   u9r  dr 

Q   =  — •  (5.23) 

I  J,     ur  dr 


in  the  physical  plane.   This  equation  can  be  used  directly 
to  evaluate  the  mean  temperature  in  domains  B  and  C,  figure 
3,  where  the  grid  lines  in  both  the  physical  and 
computational  planes  are  orthogonal. 

In  domain  A,  the  grid  lines  in  the  physical  plane  are 
oblique;  to  use  Eq.  (5.23)  will  require  a  cumbersome 
interpolation,  which  is  inconvenient.   For  this  reason,  an 
alternative  expression  is  used  as  follows: 


V  =  %1  ^  p|/z^4f  ^  d-  (5.24) 


This  expression  was  derived  based  on  an  energy  balance  of 
the  control  volume  shown  in  figure  8.  Since  only  the  inner 
wall  is  heated  in  the  present  exchanger,  the  r  value  in  the 
integrand  is  taken  to  be  unity  when  0jjj2  is  evaluated  in  the 
inner  pipe.  In  the  annulus,  this  value  is  taken  to  be  1.01 
because  of  the  dimensions  of  the  inner  pipe.  Also  notice 
that,  in  the  computational  plane,  Eq.  (5.24)  reduces  to 
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pA2^2S^ni2 


*  2  * 

pA^u^=pA2U2=pTT(r.)  u^ 


Figure  8   A  Control  Volume  Used  to  Find  the  Mean 
Temperature  by  Energy  Balance 
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The  scheme  developed  above  works  well  in  the  pipe 
regions  where  the  flow  sections  are  uniform.   At  the 
turning  tip  (point  h,  figure  3),  no  hydraulic  diameter  is 
defined.   The  characteristic  length,  C^,  at  this  point  is 
thus  arbitrarily  taken  to  be  twice  of  that  length  from  the 
tip  to  the  cap  wall,  that  is,  the  dashed  line  in  figure  2. 

Once  the  local  Nusselt  number  is  found,  the  mean 
Nusselt  number  may  be  computed  using 


/  Nu   dA    /  Nu^  r  dz 

Nu  =         ;    I, =  r-^ (5.26) 

ra     J  dA       J  r  dz 


Further  notice  that,  at  the  inlet  section,  the  local 
Nusselt  number  is  infinity  because  of  the  step  change  in 
temperature.   In  order  to  circumvent  this  problem  while 
still  evaluating  a  mean  Nusselt  number  useful  for  the 
present  investigation,  the  local  Nusselt  number  at  the 
inlet  is  estimated  by  extrapolating  those  numbers  at  two 
adjacent  points,  i=2  and  3.   The  results  are  provided  in 
the  next  chapter. 


CHAPTER  VI 
RESULTS  AND  DISCUSSION 


Numerical  results  for  the  velocity  field,  temperature 
field,  wall  vorticity,  pressure,  and  Nusselt  number  will  be 
presented  in  this  chapter  for  Reynolds  numbers  of  100,  500, 
and  1000  and  Prandtl  numbers  of  0.7,  20,  and  100. 
Computation  was  performed  on  a  Harris  800  computer,  and 
machine  plots  were  executed  with  the  use  of  a  Gould  plot- 
routine  package. 

VI. 1   Velocity  Field 
Velocity  fields  for  three  Reynolds  numbers  are 
presented  in  figure  9-   Enlarged  views  of  the  velocity 
fields  near  the  turning  point  of  the  inner  pipe  are 
provided  in  figures  10  through  12.   In  making  these  plots, 
all  arrow  heads  of  the  velocity  vectors  are  positioned 
along  lines  that  correspond  to  lines  of  constant  g    in  the 
computational  plane.   Hence,  these  arrow  heads  do  not  line 
up  radially  when  they  move  closer  to  the  turning  point 
because  they  are  in  domain  A  (figure  3)  and  in  this  domain 
these  constant  E,    lines  are  oblique  in  the  physical  plane, 
also  see  figure  4.   Certainly,  the  lengths  of  the  arrows 
still  represent  the  magnitude  of  the  velocity.   It  is  also 
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noted  that,  in  making  these  plots,  the  magnitude  of  the 
velocity  has  been  normalized  by  the  maximum  velocity  found 
for  each  Reynolds  number  tesifed.   Inasmuch  as  the  maximum 
velocity  for  Re=1000  is  larger,  the  normalized  velocities 
in  this  plot  are  smaller. 

A  close  examination  of  the  computer  data  shows  that, 
immediately  following  the  flow  entering  the  inner  pipe, 
there  is  a  slight  bulge  of  the  velocity  near  the  inner  pipe 
wall  (too  small  to  be  seen  clearly  in  figure  9).   The  fluid 
velocity  near  the  centerline,  however,  remains  uniform  at 
its  inlet  value.   Such  an  observation  appears  to  be  in 
contradiction  to  the  developing  flow  analyses  in  which  such 
a  bulge  is  not  reported.   This  inconsistency  may  be 
ascribed  to  the  fact  that,  in  the  present  investigation, 
the  axial  shear  term  in  the  momentum  equation  is  retained 
in  the  computation,  whereas  in  the  developing  flow 
analyses,  references  [18]  and  [19],  such  axial  shear  terras 
are  dropped.   Physically,  this  bulge  of  velocity  results 
from  the  need  for  flow  to  satisfy  the  continuity  equation 
and  the  no-slip  condition  simultaneously.   Because  of  the 
closeness  of  this  section  to  the  inlet,  the  momentum 
transfer  has  yet  to  reach  that  far  to  the  center  region, 
and  as  a  result,  only  velocity  nearby  is  affected  and  a 
maximum  velocity  develops  locally.   Such  a  phenomenon  has 
also  been  discussed  in  references  [20]  and  [21  J. 
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A  vortex  is  found  in  the  annulus  right  after  the  flow 
makes  the  180-degree  turn.   As  shown  in  figures  10  through 
12,  the  size  of  the  recirculation  zone  increases  with 
Reynolds  number.   For  a  quantification  of  this  observation, 
the  sizes  of  the  vortex  measured  in  distance  from  the 
turning  point  to  the  point  of  flow  reattachment  are  listed 
in  Table  1.   With  the  recirculation  region,  fluid  moves  in 
a  clockwise  loop.   The  main  flow  detaches  from  the  wall, 
resulting  in  a  depression  of  the  heat  transfer  rate  as  will 
be  discussed  later. 

It  is  also  interesting  to  note  that,  because  of  the 
blockage  of  flow  by  the  round  cap,  the  maximum  velocity  in 
the  inner  pipe  is  not  located  on  its  centerline,  but  at  a 
location  slightly  shifted  away  from  the  center.   By  the 
same  token,  the  maximum  velocity  in  the  annulus  is  not 
located  close  to  the  inner  wall,  as  is  normally  the  case 
for  an  annular  flow  [17],  but  is  displaced  slightly  tov.ard 
the  outer  wall.   This  latter  phenomenon  is  certainly  a 
result  of  the  accelerated  main  flow  that  occurs  near  the 
recirculation  region  being  displaced  away  from  the  inner 
wall . 

VI. 2   Wall  Vorticity 
The  vorticity  along  the  inner  pipe  wall  is  plotted  as 
a  function  of  position  ( dimensionless)  in  figure  13.   The 
turning  point  of  the  inner  pipe  is  located  at  0.5;  the 


Table  1 


Vortex  Sizes  for  Re=100,  500,  and  1000 
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Reynolds 
Number 

Vortex 
Size 

100 

0.37 

500 

1.37 

1000 

2.3 
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normalized  overall  length  is  1.   Again,  the  vorticity  is 
diraensionless ,  normalized  by  the  maximum  vorticity  value 
(785.105)  found  for  the  three  cases  of  Reynolds  number 
tested. 

There  are  three  places  where  the  vorticity  shows  a 
deviation:   the  inlet  region,  the  flow  region  near  the 
turning  point,  and  the  flow  in  the  recirculation  region. 
The  slight  increase  in  vorticity  in  the  inlet  region  can  be 
attributed  to  the  growth  of  the  boundary  layer  under  a  no- 
slip  condition.   A  larger  vorticity  in  this  region  enables 
these  conditions  to  be  met.   Moving  downstream,  the 
vorticity  diffuses  via  random  molecular  motion  and 
convection;  the  vorticity  tends  to  flatten  out  as  shown  in 
the  figure. 

The  turning  point  provides  a  different  perspective. 
Here  from  mathematics,  this  point  is  a  point  of 
singularity,  and  the  vorticity  approaches  infinity. 
Physically,  the  flow  accelerates  in  the  180-degree  turn, 
and  correspondingly,  the  pressure  drops  as  will  be 
discussed  later. 

As  shown  in  figure  13>  the  vorticity  changes  sign  in 
the  recirculation  region,  also  see  Table  2.   Finally,  as 
the  fluid  moves  further  downstream  toward  exit,  three 
vorticity  curves  all  merge  into  one. 
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Table  2 
Inner-Wall  Vorticity  Data  Near  the  Turning  Point 


iN.   ^^  -^ 

100 

500 

1000 

92 

15.2166 

40.7738 

62.8808 

93 

19.4822 

54.5176 

84.0731 

94 

25.5662 

74.6961 

114.966 

95 

34.3638 

105.127 

161 .401 

96 

47.5037 

152.317 

233.401 

97 

69.3520 

227.864 

346.933 

98 

114.130 

351.229 

513.606 

99 

210.448 

572.942 

786.105* 

100 

362.966 

639.621 

711  .617 

101** 

278.764 

244.890 

180.865 

102 

152.656 

1.12102 

-47.3047 

103 

26.7717 

-66.1383 

-73.1212 

104 

-0.526885 

-43.7028 

-42.0327 

105 

-6.66565 

-29.1544 

-26.0498 

106 

-7.3296 

-20.5753 

-17.5287 

107 

-6.6902 

-15.4656 

-12.7805 

108 

-5.88388 

-12.2948 

-9.94225 

109 

-5.13125 

-10.2201 

-8.10281 

110 

-4.48652 

-8,83588 

-6.86040 

*   Maximum  value 
**  Turning  point 
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VI. 3   Pressure  on  Inner  Wall 

The  pressure  distribution  on  the  inner  wall  is  plotted 
for  the  three  Reynolds  numbers  in  figure  14.   Here  the 
ordinate  is  a  dimensionless  pressure  defined  as  Eq.  (5.9); 
the  X  axis  is  again  dimensionless  as  before.   Three 
additional  informations  are  included  in  the  figure  for 
contrast.   Two  of  them  correspond  to  the  pressure  drop  for 
fully  developed  flow  in  a  circular  pipe  and  in  an  annulus, 
in  which  the  pressure  gradient  is  linear  and  is  independent 
of  the  Reynolds  number;  see  solid  lines  for  Poiseuille  flow 
and  annular  flow.   This  is  certainly  a  result  of  the 
definition  p"^ .   Another  set  of  results,  which  are  plotted 
using  the  stars  in  the  figure,  is  taken  from  the  numerical 
work  of  Hornbeck  [22]  who  evaluated  the  total  pressure  drop 
for  a  hydrodynamic  developing  flow  in  a  circular  pipe; 
again  the  slope  of  these  data  are  in  agreement  with  that 
calculated  in  the  present  exchanger.   It  should  be  noted 
that  the  level  of  these  curves  is  arbitrary;  only  the 
pressure  gradients  are  of  concern.   From  the  figure,  it  is 
clear  that  the  pressure  drop  in  the  inner  pipe  in  the 
present  exchanger  is  only  affected  near  the  turning 
point.   To  assist  this  observation,  a  dashed  line  is  drawn 
in  the  figure  to  indicate  the  approximate  location  where 
the  pressure  is  affected  in  the  inner-pipe  region. 

For  the  heat  exchanger,  there  are  three  places  where 
the  pressure  drop  appears  unusual.   Right  after  the  fluid 
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enters  the  inner  pipe,  the  pressure  gradient  is  positive. 
This  is  not  predicted  in  Hornbeck's  results  because,  in  his 
analysis,  the  axial  boundary-layer  equation  was  used  and 
the  momentum  equation  was  linearized  locally  at  any  cross 
section  z^    by  means  of  the  velocity  at  z=z-|+a2-   However, 
in  the  present  investigation,  no  linearization  is  made  and 
the  axial  diffusion  terras  are  retained  in  the  momentum 
equations.   As  a  result,  this  positive  pressure  gradient 
near  the  inlet  region  may  be  attributed  to  the  inclusion  of 
the  axial  diffusion  terms  in  these  equations.   The 
requirement  for  satisfying  the  uniform  flow  boundary 
condition  at  inlet  -Iso  contributes  to  this  pressure 
rise.   Notice  that  this  pressure  abnormality  was  also 
reported  by  Morihara  and  Cheng  [20]  for  flow  in  a  straight 
channel,  and  they  used  a  similar  mathematical  model  as  used 
in  the  present  study.   A  point  of  interest  is  that  both  the 
pressure  rise  and  the  distance  over  which  it  occurs  are 
larger  for  larger  Reynolds  numbers,  which  is  not 
unexpected. 

Before  the  flow  reaches  the  turning  point  on  the  inner 
wall,  there  is  a  rapid  drop  in  pressure.   Obviously,  this 
drop  results  from  acceleration  occurring  before  the  180- 
degree  turn  (see  dashed  line).   There  is  only  partial 
recovery  of  this  pressure  as  the  fluid  enters  the  annulus; 
as  would  be  expected,  the  pressure  loss  is  partially 
irreversible . 
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The  circulation  of  flow  in  the  annulus  gives  rise  to 
another  region  of  positive  pressure  gradient.   Finally,  all 
three  curves  merge  toward  the  exit  of  the  exchanger  in 
satisfaction  of  the  conditions  imposed  there,  Eqs.  (5.4) 
and  (5.13).   The  flow  becomes  fully  developed  as  evidenced 
by  the  solid  line  drawn  in  the  figure.   To  test  the 
accuracy  of  the  pressure  computed  in  the  numerical  work, 
dpgj,/dz  is  evaluated  at  the  exit  section.   As  listed  in 
Table  3>  this  pressure  gradient  has  a  value  ranging  from 
-1.0245  to  -1.0781;  the  error  is  thus  about  8%   as  compared 
with  the  pressure  gradient  for  annular  in  a  fully  developed 
condition;  also  see  Eq.  (5.13). 

VI .4   Temperature  Field 
Isotherms  are  plotted  in  figures  15  through  17. 
Recall  that,  for  the  problems  under  investigation,  0=0 
along  the  inner  wall,  and  9^=0  along  the  outer  wall. 
Notice  that,  each  figure  presents  results  for  all  three 
Reynolds  numbers  at  one  Prandtl  number.   Hence,  comparison 
between  the  isotherm  distributions  of  one  figure  provides 
insight  into  the  effect  of  Reynolds  number;  a  cross 
reference  of  the  plots  at  the  same  Reynolds  number  offers 
an  assessment  of  the  effect  of  the  Prandtl  number. 

The  boundary  conditions  required  that  all  isotherms 
emerge  from  the  junction  of  the  inner  wall  with  the 
inlet.   When  one  compares  any  position  downstream  of  the 


Table  3 
Dlniensionless  Pressure  Gradients  at  Exit 


68 


Re 

+      * 
dp   /dz 

ex 

100 

-1  .0245 

500 

-1 .0781 

1000 

-1.0687 

*  Correct  value:   -1 
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entrance,  it  is  observed  that  an  increase  in  Reynolds 
number  at  the  same  Prandtl  number  produces  a  steeper 
temperature  gradient  near  the*  (inner)  wall. 

A  cross  reference  of  the  plots  in  the  three  figures  at 
fixed  Reynolds  number  reveals  that  an  increase  in  the 
Prandtl  number  produces  the  same  general  effect  as  an 
increase  in  Reynolds  numbers.   This  is  to  be  expected  since 
increased  Prandtl  number  implicates  a  smaller  thermal 
diffusion  as  compared  to  the  momentum  diffusion. 
Therefore,  heat  entering  the  fluid  at  the  inner  wall  is 
unable  to  penetrate  as  deeply  into  the  fluid.   In  fact, 
figure  16  (a)  shows  these  isotherms  are  so  crowded  near  the 
wall  that  full  presentation  of  them  becomes  impossible; 
only  those  isotherms  from  0.9  to  0.99999  are  thus  shown. 

For  all  cases  tested,  the  isotherms  are  crowded  toward 
the  turning  point  of  the  inner  pipe.   This  results  from  the 
high  velocity  in  this  region.   With  respect  to  the  large 
wall  gradients  isotherms  observed  near  the  inlet  section 
and  near  the  turning  point,  it  should  be  noted  that  these 
imply  high  convective  coefficients  in  these  regions. 

A  totally  different  state  of  affairs  occurs  in  the 
recirculation  region  following  the  sharp  turn.   Figure  15 
(b)  and  (c)  show  a  smaller  temperature  gradient  in  these 
regions,  indicating  a  lower  heat  transfer  coefficient. 
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Further  downstream  where  the  fluid  reattaches  to  the 
wall,  a  crowding  of  the  isotherms;  see  figure  15  (b)  and 
(c),  indicates  slight  rise  of  the  convective  coefficient. 

These  observations  are  summarized  in  a  series  of 
Nusselt-number  plots  presented  in  the  next  section. 

VI. 5   Local  Nusselt  Number 


Local  Nusselt  numbers  for  the  nine  cases  are  shown  in 
figures  18  through  20  at  three  Prandtl  numbers.   Two  pieces 
of  information  are  also  included  in  figure  13;  one  relates 
to  the  Nusselt  number  for  fully  developed  flow  in  a 
circular  pipe  under  a  constant-wall-temperature 
condition.   The  Nusselt  number  for  this  case  is  3-658  [18]; 
see  horizontal  line  drawn.   The  other  set  of  data  is  taken 
from  the  numerical  work  of  Hornbeck  [23]  who  studied  the 
Nusselt  number  under  a  developing  flow  and  heat  transfer 
condition  at  the  same  Prandtl  number  (0.7).   His  results 
are  plotted  using  stars  in  the  figure. 

Referring  to  figure  18,  it  is  seen  that  the  present 
Nu^  values  are  in  good  agreement  with  the  results  of 
Hornbeck  for  Re=100.   Moving  downstream,  the  local  Nusselt 
numbers  asymptotically  approach  those  for  fully  developed 
flow.   Near  the  turning  point,  the  turn  takes  effect, 
resulting  in  an  increase  in  heat  transfer  and  the  local 
Nusselt  numbers  increase  rapidly.   Again,  to  assist  this 
observation,  a  dashed  line  is  drawn  in  this  figure  to 
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indicate  the  approximate  location  where  turn  takes 
effect.   For  the  higher  Reynolds  number  case,  the  local 
Nusselt  numbers  near  the  inlet  region  of  the  present  study- 
deviate  only  slightly  from  the  Hornbeck's  results. 

For  all  the  cases  studied  for  the  heat  exchanger,  a 
large  Nusselt  number  is  found  near  the  inlet,  the  turning 
point,  and  in  the  reattachment  region  and  low  Nusselt 
number  occurs  in  the  recirculation  region. 

It  is  clear  that  once  the  flow  is  deflected  into  the 
annulus,  a  redevelopment  of  flow  occurs.   In  the  annulus 
region,  a  fully  developed  condition  is  reached  earlier  for 
fluid  of  low  Peclet  number,  see  figures  18  through  20.   In 
order  to  check  on  the  computations  of  the  present  work, 
Lundberg  and  coauthors'  results  [24]  are  used.   As  shown  in 
Table  4,  Lundberg's  results  for  fully  developed  velocity 
and  temperature  profiles  in  an  annulus  appear  to  be  in  good 
agreement  with  the  present  investigation  that  is  computed 
at  Pr=0.7  and  Re=100.   If  a  linear  interpolation  based  on 
Lundberg's  results  at  r"*'=0.5  and  1.0  can  serve  as  a  guide, 
then  the  error  is  only  1.5^. 

VI. 6  Mean  Nusselt  Number 


Figure  21  presents  mean  Nusselt  number  as  a  function 
of  Reynolds  number  with  Prandtl  number  as  a  parameter. 
This  figure  shows  that  mean  Nusselt  number  increases  with 
both  the  Reynolds  and  Prandtl  numbers.   A  correlation  of 
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Table  4 

Comparison  of  the  Computed  Local  Nusselt  Number 
at  the  Exit  with  Lundberg's  Data 


+ 
r 

Nusselt  Number 

0.5* 

5.738 

0.711** 

5.297 

1.0* 

4.861 

*      Lundberg's  data 

**  Present  investigation  for  Pr=0.7  and 
Re=100 
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10^  I- 
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3x10 


3x10^ 


Re 


Figure  21   A  Plot  of  Mean  Nusselt  Number 
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the  results  yields 

Nu  =  1.15RB°-528pr0-2S8  (6.1) 


m 


which  is  valid  for  Re=100  to  1000,  and  Pr=0.7  to  100. 
Variance  evaluated  based  on  (n-2)  degree  of  freedom  (n 
represents  number  of  correlated  point)  is  2.6  [25]. 

The  present  heat  exchanger  is  found  to  be  far  more 
effective  in  heat  transfer.   To  verify  this  point, 
Hornbeck's  data  [23]  for  air  (Pr=0.7)  under  a  developing 
velocity  and  temperature  condition  in  a  circular  pipe  are 
used  for  comparison.   As  listed  in  Table  5,  the  mean 
Nusselt  number  for  the  SPRF  heat  exchanger  is  1.57  times 
that  for  pipe  flow  when  the  Reynolds  number  is  500. 
Doubling  this  Reynolds  number  enables  this  Nusselt-number 
ratio  to  be  raised  to  1.87.   The  superiority  of  the  present 
exchanger  is  thus  clearly  seen. 

VI. 7   Truncation  Error  and  Convergence 
In  order  to  check  the  truncation  error  and  the 
convergence  of  the  solutions  presented,  three  different 
grid  systems  are  used.   It  is  noted  that,  for  the  ease  of 
comparison,  the  number  of  grid  points  along  the  constant  r\- 
line  is  fixed  at  163;  the  number  of  grid  points  along  the 
constant  E,-line    is  tested  for  21,  41,  and  61.   This  results 
in  lines  of  constant  c  for  these  grid  systems  having  the 
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Table  5 

Comparison  of  Mean  Nusselt  Number  Between 
SPRF  Heat  Exchanger  and  Heating  in  a  Circular  Pipe 


Pr 

Re 

Mean  Nusselt  Number 

Mean  Nusselt- 
Number  Ratio 

Hornbeck's  Data* 

Present  Study 

0.7 

500 

7.4 

11.62 

1.57 

0.7 

1000 

9.57 

17.88 

1  .87 

*  Developing  velocity  and  temperature  profiles  in  a 
circular  pipe 
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same  contours  in  the  physical  plane  but  with  different  grid 
sizes . 

Figure  22  shows  the  comprtited  stream  functions  along 
constant  5-lines  for  the  three  grid  systems.   In  these 
figures,  i=50  and  152  correspond  to  the  mid-section  in  the 
inner  pipe  and  the  annular  region,  respectively;  i=101 
corresponds  to  the  constant  ^-line  at  the  turning  point 
(see  figure  5).   It  is  clear  that  the  computed  stream 
function  values  based  on  the  use  of  grids  163x41  converge 
to  the  results  of  the  dense  grid  (153x61).   A  close 
examination  of  the  difference  between  these  two  sets  of 
data  shows  that  the  maximum  deviation  is  only  about  3% 
relative  to  the  dense  grid.   However,  the  CPU  time  required 
for  computation  with  the  dense  grid  is  increased  by  a 
factor  of  about  3.5  (Table  6)  using  the  Harris  800  computer 
on  campus.   For  this  reason,  a  denser  grid  is  not  used  for 
the  present  computation.   The  results  for  163x41  are 
considered  converged. 
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63x61) 
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inside  the  annular  region) 
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0.0 


(163x21) 


0. 


(163x41) 


3   .4    .5   .6    .7   .8 
(^-^min^/^Vx-^nin) 


Figure  22   Comparison  of  the  Computed  Stream  Functions 

Along  Constant  C-line  for  Different  Grid  Systems 
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Table  6 

CPU  Time  Required  for  Computation  of 
Stream  Functions  Given  in  Figure  22 


\^^--^^^   Re 

\  CPU     *""^^ 

\   time 

Grids  \ 

100 

500 

1000 

163x21 

1  hr 
30  rain 

57  rain 

1  hr 
17  rain 

163x41 

4  hr 
55  rain 

4  hr 
1 4  rain 

6  hr 
19  rain 

163x61 

12  hr 
15  rain 

10  hr 
3  rain 

18  hr 
16  min 

*  CPU  time  is  based  on  the  use  of  Harris  800 
computer  on  campus 


CHAPTER  VII 
CONCLUSIONS  AND  RECOMMENDATIONS 


A  numerical  solution  is  given  for  studying  the  flow 
and  heat  transfer  in  a  single-pass,  return-flow  heat 
exchanger.   Stream  function  and  vorticity  are  used  in  the 
formulation,  and  the  original  irregular  geometry  of  the 
exchanger  in  the  physical  plane  is  transformed  into  a 
rectangular  domain  with  square  grids  in  a  computational 
plane.   Elliptic  partial  differential  equations  are  used  in 
the  transformation,  and  grid  lines  are  clustered  with  the 
use  of  separate  clustering  functions  which  allow  fine  grids 
to  be  formed  in  regions  where  the  fiov7  and  temperature 
fields  change  most  rapidly.   The  source  terms  in  the 
Poisson's  equations  are  also  evaluated  in  the  solution. 

A  finite  difference  method  is  used  to  solve  the 
problem  in  the  computational  plane.   In  the  numerical 
algorithm,  in  order  to  assure  a  convergent  solution  in  the 
iteration,  the  diagonal  terms  in  the  resulting  matrix 
equations  are  reconstructed  to  make  them  dominant. 
Vorticity  boundary  conditions  at  the  walls  are  expressed  in 
a  first-order  approximation.   The  velocity  and  temperature 
conditions  along  the  centerline  of  the  heat  exchanger  are 
also  finite  differenced,  with  the  first-order  differential 
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terms  of  velocity  and  temperature  expressed  in  a  first- 
order  approximation.   Solutions  are  obtained  for  Pr=0.7, 
20,  and  100,  and  Re=100,  500,*  and  1000. 

Numerical  data  show  a  marked  increase  in  velocity  as 
flow  makes  a  180-degree  turn  in  the  end  cap.   Here  a  rapid 
drop  of  pressure  is  also  noticed.   Once  the  flow  enters 
annulus,  the  main  flow  is  separated  from  the  inner  wall, 
and  a  recirculation  is  found  near  the  wall.   This  results 
in  a  drop  in  the  heat  transfer  coefficient  in  this 
recirculation  region.   Moving  further  downstream,  the  flow 
reattaches  to  the  wall,  causing  a  slight  rise  in  the  heat 
transfer.   The  mean  heat  transfer  coefficient  for  the  heat 
exchanger  is  found  to  be  proportional  to  the  Reynolds 
number  and  Prandtl  number.   The  numerical  results  are 
correlated  to  develop  an  empirical  equation. 

For  a  check  of  the  accuracy  in  the  solution,  the 
pressure  gradient  and  the  Nusselt  number  at  the  exit  of  the 
exchanger  are  also  evaluated  and  compared  with  the  theory 
(analytical  equation)  and  practice  (literature  data).   The 
agreement  is  satisfactory.   The  single-pass,  return-flow 
heat  exchanger  is  found  to  have  a  mean  Nusselt  number  that 
is  about  1.57  to  1.37  times  that  for  heat  transfer  in  a 
singular  pipe  under  comparable  working  conditions  (Pr=0.7, 
and  Re=500,  and  1000).   The  exchanger  is  thus  shown  to  be 
effective  in  heat  transfer. 
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Based  on  the  study  made  in  this  dissertation, 
directions  of  further  research  are  charted  as  follows. 

1.  There  are  a  wide  range  of  variations  in  the 
conditions  used  which  the  SPRF  exchanger  can  be 
operated.   They  include  (i)   change  in  the  flow 
direction,  (ii)   change  in  the  system  geometry, 
particularly  that  gap  size  from  the  inner  pipe  to 
the  end  cap,  (iii)   change  in  thermal  boundary 
conditions  on  both  walls,  and  (iv)   introducing 
minor  modifications  in  the  construction  of  the 
exchanger,  such  as  eccentricity  in  the  pipe 
layout,  flow  perturbation  devices,  etc.   They 
provide  a  wide  range  of  conditions  to  be  tested  in 
the  future. 

2.  Because  of  the  complexity  of  the  geometry  of  the 
exchanger,  it  provides  an  opportunity  for  testing 
various  numerical  methods,  particularly  those 
approximation  methods  in  the  finite-difference 
formulation  of  the  boundary  conditions. 

3.  A  study  is  not  complete  without  experiment.   It  is 
highly  desirable  that  hardware  be  built  to  test 
the  results  obtained  in  the  numerical  solution.   A 
limited  experimental  program  serves  well  to 
validate  the  numerical  solution;  until  then  the 
numerical  solution  can  be  used  with  confidence  for 
acquiring  data  for  various  operating  conditions. 


APPENDIX  A 

DERIVATION  OF  TRANSFORMATION  RELATIONS, 

GOVERNING  EQUATIONS,  AND  BOUNDARY  CONDITIONS 

IN  THE  TRANSFORMED  PLANE 


This  appendix  gives  the  derivation  of  the 
transformation  relations,  governing  equations,  and  boundary 
conditions  in  the  physical  and  computational  planes  when 
elliptic  partial  differential  equations  are  used  to 
generate  the  computational  grid  system. 

Referring  to  figure  A.1,  the  independent  variables  in 
the  physical  plane  are  z(C,n)  and  r(C,n),  while  those  in 
the  computational  plane  are  ?(z,r)  and  n(z,r).   In  the 
computational  plane,  the  domain  of  interest  is  rectangular 
in  shape  and  grids  in  this  plane  are  uniformly  spaced 
(i.e.,  square  grid).   Since  all  the  computations  are  to  be 
made  in  the  transformed  domain,  the  independent  variables 
in  the  governing  equations  and  boundary  conditions  must  be 
changed  from  the  physical  domain  to  the  transformed  domain. 

Consider  the  vector 


z(5,n) 
r(5,n) 


(A.1) 


which  may  be  differentiated  to  give 
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(A. 2) 


Cramer's  rule  may  be  used  to  solve  for  d?  and  dn  as 


d? 
dn 


.  n 
— r 


5    5 


dz 
dr 


(A. 3) 


where  J  is  the  transformation  Jacobian 


J  = 


Similarly, 


9(2, r)  ^ 
3C5,n) 


z_  z 

C  n 

r  _  r 

5  n 


dS 
dn 


?z   ?r 
^z   ^r 


dz 
dr 


(A. 4) 


The  following  identities  may  be  established  by  comparing 
(A. 3)  and  (A. 4) 


C   =  r  /J 
z    n 


5r  -  -  ^« 


"z  =  -  'i'' 


(A. 5) 


and 


^r  =  \/^ 
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The  original  partial  differentials  in  the  physical 
plane  may  be  transformed  to  the  computational  plane  by 
using  the  chain  rule  for  partial  derivatives.   It  follows 
that 


dz         ^z    35    ^z  8ti 

=  (r  ^  -  r  ^)/J  (A. 6) 


3r    ^r  95    ^r  3n 

=  (-Z  -^  +  z^  -^)/J  (A. 7) 

Ti  35     5  9n 


2 

1 =  ^-  (r   -2- 

^^2    3Z  ^z  35  '  "z  3ti' 

_      _3 _3 /•_3 \        _3 _3 /_3 s 

^ZZ  35    ^z  3Z  35     '^zz  3ri    ^z  3Z  3n 

=  5    ^-  +  5  [  (  5   ^-  +  n   -^)-^] 
'^zz  35    ^z^'^z  35    ^z  3n  8C 

o   2  „   2 

^zz  35    '^z        2        ^zz  3n    ^z    2 
35  3ri 

+  2  5  n  ^ (A. 8) 

^z"z  358n 


and 
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3^     ^    9      ?  8^         a      2  8^ 


r  r  95  on 


Combining  (A. 8)  and  (A. 9)  gives  the  Laplace  operator 


^^     ^^   =  (C?  +  C?)  ^  +  2(?  n^  +  £  n  )  ^ 


„  2  ^  2    '  z     r'  .^2     '"z  z    "r  r'  a^Sn 
9z    or  35 

^2    2s  3^     /  ^      r   ^  3 

^  ^^z  ^^^7-2"  ^^zz  "  ^rr^  T? 
3ri 

+  (n    +  n   )  4-  (A. 10) 

zz    rr   3n 


in  which  the  coefficients  for  partial  differentials  on  the 
right-hand  side  may  be  expressed  in  terms  of  J  as 


and 


?2  +  e^  =  (z^  +  r2)/j2  =  a/j2 
z    r    ^  n    n 


2        2 

C  n   +  5  n   =  -  (r^r   +  z-z  )/J   =  -  b/J 
z  z    r  r       C  n    5  n 


2    2,2    2s/t2     /,2 
ri„  +  n„  =  (z^  +  r  )/J   =  Y/J 


where 


2    2 

a  =  z   +  r 
n    n 


=  r^r^  +  z_z 


n    ?  n 


2  2 


Then, 
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2             2                         2 
9         ^  i 1     ^      J o„    9 

+   ?r  =   -7-   (>  a  5 2B  — 


.2         ,    2    ~    J 
9z  3r 


35 


953n  9ti 


^  ^zz         ^rr'    dE,  'zz         ^rr^    an 


(A. 11 ) 


If  Laplace  equations 


C    +5    =0 
zz     rr 


(A. 12) 


and 


0 


(A. 13) 


are    used    in    the    transformation,    Eq.    (A.11)    may   be 
simplified    to   be 


1/8  ,„    af ^ 


2    ■  2         .2    '"         2 

3Z  9r  J  95 


Y  — y. 
3  5  9ti  9ri 


(A. 14) 


If  Poisson  equations 


?zz  ^  S^r  =  P^^'^^ 


(A. 15) 


and 


^zz  ^  ^rr  =  ^^^'^) 


(A. 16) 


are  used,  Eq.  (A.11)  reduces  to 


94 


2  2  2  2  2 

\  a  —  -    2  e   +    Y   ^) 


p  P    ~      P  2  ? 

3z  3r  J  3C  859n  3ti 

+  P  TT  +   Q  T^  (A. 17) 


Those  relations  derived  above  can  be  used  to  derive 
transformation  equations  listed  below. 

(1)   From  Eq.  (A. 14),  the  Laplace  transformation 
equations  in  the  computational  plane  become 


az^^  -  2B2^   +  yz    =0  (A. 18) 

C?      5n     nn 

and 

ar^^  -  26r^   +  yr    =0  (A. 19) 

£,i  5ri      nn 


(2)   From  Eq.  (A. 17),  the  Poisson  transformation 
equations  in  the  computational  plane  become 


az    -  2bz    +  vZ    =  -  J^(Pz   +  Qz  )  (A. 20) 

5C      5n     nn         5     n 

and 

ar    -  2Br    +  ^r    =  -  J^(Pr^  +  Qr  )  (A. 21) 

CS      5n     nn  5     n 

From  these  two  equations,  the  source  terms  P  and  Q  can  be 
derived  as 


P  =  (z  Dr  -  r  Dz)/J^  (A. 22) 

n      n 

and 
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Q  =  (r^Dz  -  z^Dr)/J^  (A. 23) 

where  Dr  =  ar^^  -  26r^^  +  Yr^^  and  Dz  =  az^^  -  26z^^  +  Yz^^, 

(3)   With  the  use  of  Eqs.  (A. 7)  and  (A. 17),  the 
vorticity  definition  equation  (2.2)  in  the  computational 
plane  takes  the  following  form 


Jz  Jz 

=  -J^r^  (A. 24) 


(4)   With  the  use  of  Eqs.  (A. 6)  and  (A. 7),  the 
velocity  equations  (2.3)  and  (2.4)  in  the  computational 
plane  become 


^  =  jF  ^^i\    -   ^nh^  ^^-25) 


and 


^  -JF  ^H\  -   ^^)  ^^-26) 


(5)   With  the  use  of  Eqs.  (A. 6),  (A. 7),  and  (A. 17), 
the  vorticity  transport  equation  (2.1)  in  the  computational 
plane  becomes 


2  Jz  Jz 

aa^^   -  26fi^^  +  yn^^  -  ^  a  ^  (j2p  _  -^)^^   +  (j^q  +  -^)^^ 

r 

=  ^  [i\^,    -    ^,\)    ^f  (r,>^5  -  r^^^)]  (A. 27) 
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(6)   The  energy  equation  (2.5)  in  the  computational 
plane  is  also  derived  as 

«^5  -  230^,  .  ye^,  .  (J2p  -  ^)e^  .  (j^q  .  ^)e^ 

Derivation  of  the  governing  equations  in  the 
computational  plane  is  now  complete.   To  derive  boundary 
conditions,  one  must  first  derive  the  expressions  for  the 
normal  and  tangential  vectors  in  the  computational  plane. 

For  a  scalar  function  F  in  the  physical  plane,  the 
gradient  of  F  is 


VF  =  F  e   +  F  e  (A. 29) 

z  z    r  r 


This  equation  can  be  recast,  with  the  use  of  (A. 6)  and 
(A. 7),  as 


VF 


=•  7  ^^\h   -   'i'/'z   *    (-/?  ^  ^/n>«ri     (*-30) 


On  a  line  of  constant  n,  n^  =  0  and  n^  =  1 ;  Eq.  (A. 30) 
reduces  to 


„  .  1.  (-,^e^  .  z^e^)  (A.31) 
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Similarly,  on  a  line  of  constant  5,  5^  =1,  5^  =0,  and 


'^   =T^\K  -  \K^  (A-32) 


With  reference  to  figure  A. 2,  the  unit  normal  vector 
to  lines  of  constant  n  can  be  expressed  as 

n^^)  -^-        ^  -    ^  -  (A.33) 

The  unit  normal  vector  to  lines  of  constant  5  is 


^  /    \  r  e   -  z  e 


The  unit  tangent  vectors  can  be  derived  with  the  use  of 
cross  products  as  follows: 


(n)   :(n)v  :   _  ^g^z  ^  ^K^r 


t^^^  =  n^^^X  e,  =  -^-^ ^-^  (A. 35) 


and 


The  directional  derivatives  of  F  along  normal  and 
tangential  directions  can  be  expressed  as 


~^^ms^' 


<    QJ 


-P 

c 

CO 

•p 

CO 

c 
o 
u 

u 
o 

CQ 


c 
o 

m 
u 
o 
-p 
o 
<u 
> 


CO 

s 
ti 
o 

3 

T3 
C 
CO 

-P 
C 

<U 
bO 
C 
CO 

EH 


CM 
< 
0) 

u 

:3 

bO 


99 


8n 


9F     ;(Ti) 


3F   _  -(n) 


:i^ 


8F     ^(Z 
-   n 


7T? 

on 


and 


9F    ;(?) 


VF  =  4 ^^ 


r^ 


VF 


rr 


r^ 


VF 


(A. 37) 


(A. 38) 


(A. 39) 


(A. 40) 


The  arc  length  dt  along  any  contour  r  in  the  (z,r^ 
plane  is  (see  figure  A. 3) 


dt 


(r 


^  =  \    (dz)2  +  (dr)2 


In  the  transformed  plane, 


dt 


(r) 


-  J^(dC) 


^  +  BdCdn  +  a(dn)^ 


(A. 41) 


(A. 42) 


where  (A. 2)  has  been  used  for  dz  and  dr. 

If  r  is  a  line  of  constant  n,  Eq.  (A. 42)  reduces  to 


,(n) 


dt'  '  =  \^(      d? 
Similarly,  if  r  is  a  line  of  constant  C,  then 


(A. 43) 
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Figure  A. 3   Schematic  Showing  the  Arc  Length  Along 
Contour  r 
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Boundary  conditions  will  now  be  derived  in  the 

transformed  plane.   Note  that  all  conditions  of  the  first 

kind  can  be  used  as  is;  no  transformation  is  necessary  of 

them.   As  for  the  other  conditions,  notice  the  fact  that, 

in  both  inlet  and  exit  domains  B  and  C,  grids  are 

orthogonal  in  both  original  and  transformed  coordinates. 

So,  metrics  such  as  r_,  z  ,  r^..,  and  z   all  vanish. 

Furthermore,  terms  such  as  r^   and  r  ^  also  vanish  because 

5n      n5 


r    =  -^(-2^)  =  -^(^^)  =  r  (A. 45) 


and  r   =0. 

With  these  in  mind,  consider  Eq.  (2.7a)  first.   It  is 
clear  that  this  condition  is  equivalent  to 

u   =v   =;h   =n   =0  (A. 45) 

Next  consider  Eq.  (2.7b).   To  express  it  in  the  trans- 
formed plane  requires  use  of  Eq.  (A. 8)  in  which  the  last 
three  partial  differentials  all  vanish  because  n^  =  0*   '^^^ 
remaining  two  terms  on  the  right-hand  side  of  this  equation 
may  be  transformed  by  use  of  Eqs.  (A. 5)  and  (A. 6).   It 
follows  that 
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Q    =^0^  (A. 47) 


Along  the  axis  of  symmetry,  that  condition  Uj,  =  0  is 
changed  to 


z  u   -  z  u^  =  0  (A. 48) 


using  (A. 7).   Similarly,  e   =  0  becomes 


z,e  -  z  e^  =  0  (A. 49) 


At  the  outer  pipe  wall,  e^  =  0  is  changed  to 


y9  -  ee,  =  0  (A. 50) 


because  of  (A. 37). 

The  vorticity  boundary  conditions  at  the  inner  and 
outer  pipe  walls  can  be  derived  using  a  lengthy  procedure 
as  follows. 

At  the  wall,  u  =  v  =  0;  Eqs.  (4-3)  and  (4.4)  become 


z    i>      -  z  J;   =0  (A. 51  ) 


and 


'^\    -  r^*^  =  0  (A. 52) 


These  walls  may  be  regarded  as  streamlines  along  which 
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4-^=0  (A. 53) 


everywhere . 

Using  this  relation  in  (A. 51)  yields 


l*^  =  0  (A. 54) 


which  is  also  valid  everywhere  along  the  wall. 

Differentiating  Eq.  (4.3)  with  respect  to  5  gives 


(>A.55J 


(Jr)2 


where  the  second  term  in  the  numerator  vanishes  because  of 
(A. 53)  and  (A. 54).   It  follows  that 


u      ^  -r-   (z  ^   +  z  ij;    -  z   i>   -  z J>      )  (A. 56) 


Along  the  wall,  u^.  =  0.   Also  notice  that  (A. 53)  can 
be  used  to  write 


^l>^^    =   0  (A. 57) 


Then 


from  Eq.  (A. 56). 


4/,   =  0  (A. 58) 
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Equations  (A. 53),  (A. 54),  (A. 57),  and  (A. 58)  may  now 
be  used  in  (4.2)  to  write  ^  at  the  wall  as 


n   =  -  -I-  ^^^  (A. 59) 

w     j2^   nn 


APPENDIX  B 
DERIVATION  AND  USE  OF  EQUATION  (3.6) 


This  appendix  gives  the  details  on  how  Eq.  (3.6)  is 
derived  and  used  to  locate  grid  lines  in  z  direction. 

It  is  necessary  to  locate  the  grid  lines  in  figure 
(B.1)  so  that  Zq  corresponds  to  index  ig,  z-j  to  i-^,    Zj_  to 
ij^,Z£_^  to  if-1  ,  and  finally  z^  to  i^.   In  particular, 
because  of  the  z^  function  to  be  chosen  later,  it  is 
necessary  that  the  following  conditions  be  met 


for  i.   =  i^  (B.1) 

z^  =  z^_^      for  i^   =  i^-1  (B.2) 

z.  =  z„        for  i.   =  i„  (B.3) 
if             if 


In  order  to  meet  the  other  conditions  mentioned,  but  not 
formulated  above,  z^    is  chosen  to  be  a  cubic  polynomial  of 
the  form 


z.  =   Zq    +   a.'?.  +  a^'V^    +   a^'V^  (B.4) 

where 

i .  -  i^ 


1    1^  -  1q 
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There  are  three  unknown  coefficients  (a>|,  a2,  a^)  in 
Eq.  (B.4)  and  there  are  three  conditions  ((B.1),  (B.2), 
(B.3))  to  solve  them.   These  coefficients  are  found  to  be 


a.  =  z^  -  Zq  -  a2  -  a^  (B.5) 

a^  =  [Az^  -  h(z^  -  Zq)  -  a^(h^  -  h)]/(h^  -  h)      (B.6) 
a   =  [Az^  +  Az^  -  2h(z^  -  ZQ)]/(h  -  3h^  +  2h^)     (B.7) 


where 


AZ^     =    Z^     -    Zq 


AZ^  =  2f  -  ^t.^ 


h  =  (i^  -  Xq)' 


In  practice,  values  of  Iq,  i^,  Zq,  z-^,  Az^  and  Az^ 
must  be  preassigned.   With  suitable  values  of  Az^  and  Az^ 
(to  be  tested),  Eq.  (B.4)  yields  a  monotonous  function  of 
increasing  grid  sizes  from  Zq  to  z^. 


APPENDIX  C 
DERIVATION  OF  FINITE  DIFFERENCE  EQUATIONS 


Transformation  equations  (3.3)  and  (3.4)  and  governing 
equations  (4.1)  and  (4.5)  will  be  changed  to  finite- 
difference  forms  in  this  appendix.   For  the  ease  of 
formulation  of  these  equations  and  later  computation  in  the 
transformed  plane,  the  grids  in  the  computational  plane  are 
taken  to  be  squares,  i  .e .  ,  AC  =  An  =  1 .   In  the  indexing 
scheme,  i  refers  to  5  and  j,  n  axis. 

Defining  F  =  F(c,n)  as  a  twice  dif ferentiable  function 
of  5  and  n,  the  finite  difference  representations  for  the 
first  and  second  derivatives  of  F  at  an  interior  point 
(i,j)  may  be  expressed,  by  means  of  a  second-order  central 
difference  scheme,  as  follows: 


(^c'l.a  "<^.1, 0-^-1. a>/2  '^•^' 

(^«)i,j  =  ^1*1,0  -2^,0  *^i-i,d  (  =  -3) 

(F   )     =  (F        -  F        +  F 

^  CTi^i,j    ^  i  +  1,j  +  1     1  +  1, j-1  ^   1-1, j-1 

-  ^-1, 0  +  1^/4  (C.4) 

and 
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(^.>l,a  =  "i,J*i  -  ""i,J  *   'i,3-i 
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(C.5) 


The  first-order  backward  and  forward  differences  of  F 
and  F  are  given  respectively  as 


and 


(F  ).  .  =  F.  ,  .  -  F.  .                          (C.8) 

^  c'i,a  1+1, J   i,j 

(F  ).  .  =  F.  .  ,  -  F.  .                          (C.9) 

n  1,3  1,0+1     1,0 


To  facilitate  representation  of  partial  differentials 
in  computer  programs,  the  changes  of  F  in  both  g   and  ri 
directions  are  written  compactly  as 


F   =  F.  ,  .  -  F.  ,  .  (CIO) 

5     1+1, J     1-1,0 

F   =  F.  .  ,  -  F.  .  .  (C.11) 

and 

^^Ti  ^  ^1  +  1,0  +  1  "  ^i  +  1,0-1  ■"  ^i-1,d-1 


Then,  the  transformation  parameters  a,  3,  y,  J,  P»  3.nd 
become 
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and 


where 


and 


a  =  (i2  ,  ?2)/4 


e  =  (5^i^  +  ?g?n)/4* 


Y  =  (5^  +  ?^)/4 


J  =  (^c^n  -^  ^n^?)/^ 


P  =  (i^Dr  -  ?,Dz)/(2J^) 


=  (r^Dz  -  z^Dr)/(2J^) 


Dr  =  -iv^^).^.    -   26(r,,),^.  .  Y(r^^),^. 


Dz  =  <^(z,,)i,j  -  23(z^^),^.  .  T(z^^),^. 


The  transformation  equations  (3.3)  and  (3.4)  take  the 
following  forms  in  the  computational  plane 


az.  .  .  -  2(a+y)z.   .  +  az .     . 
1-1, J  i,J      1+1, J 


-•"^1,0*1  *  ^i,a-i*  *2^5n        'c-iJ) 


and 


cir.  .  .  -  2(ci  +  y)r.  .  +  "r.  .  . 
1-1,0  i»a     1+1, J 


-''('■l.a.l  *  'i.3-^^   *  2  ?5n  (C-^t' 
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These  equations  are  to  be  used  to  map  the  grid  points  from 
the  computational  plane  to  the  physical  plane. 

With  the  finite  difference  equations  given  above,  the 
governing  equations  in  Chapter  IV  can  be  changed  to  finite- 
difference  forms  as  follows: 

Vorticity  transport  equation: 


[2(a+T)  +  ^  +  Re(— ^  ^^ 1-  ij;  )J 


4p2   c    ^^2   n    i,j 


1  /  .2 


Jz 


1  /  .2, 


J: 


J~  Jz 


-  R^  4F  ^^5  '  ^^  4F  ^c%  "  2  ~\r. 


(C.15) 


Vorticity  definition  equation: 

-    U    .  1<J2P  .  ^)]*i,i,j  *  La  -  1(J2P  .  ^)]ti.i,j 
Jz  Jz 


+  J  rn.  .    -  ^  ^^ 
1 ) J    2  ^5n 


(C.16) 


Velocity  equations 


^,j  =  4Jr  "^n  "  4Jr  "^s 


(C.17) 
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TT^^   -  mrr^.  CC.18) 


1,0  -  4Jr  %   4Jr  ^5 

t 

Energy  equation: 
2(ci+Y)e.  . 


It  is  convenient  to  introduce  the  following 
expressions  for  metric  coefficients 

CO  =  -  b/2 
CI  =  2(a+Y) 

P     Jz 
C2  =  a  +  (J^P  +  -2^)/2 


C3  =  a  -  (J^P  +  -2^)/2 


P     Jz 
C4  =  y  +  (J^Q  -  -2f)/2 


C5  =  Y  -  (J  Q  -  -2r^/2 


C6  =  J^r 


C7  =  2(a+Y)  +  J^/r^ 
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C8    =    a    +    (J^P 2F^/2 


C9    =    a    -    (j2p    _  ^)/2 


CIO    =    Y    +    (J^Q    +  -2f)/2 


?  Jz 

C11    =   Y    -    (J'^Q   +  -2f)/2 


C12   =   z^/(4Jr) 


C13    =   i^/(4Jr) 


C14    =   r^/(4Jr) 


C15    =   r^/(4Jr) 


and 


C16   -  J/(4r) 


C17   =    (J?^)/(4r^) 


CIS   =    (J?^)/(4r^) 


so    that   Eqs.    (C.15)    to    (C.19)    can   be    recast   compactly  as 


[C7    +   Re(C17K    -   C18^    )]fi.     . 


-   ReCl6(?    n^    -    lb    n    )    +   COf^, 
n    E,         ' E,    n  Sn 


(C.20) 
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+  C6j2.  .  +  CO^^    •  (C.21) 

u    =  C12^   -  C134i,  (C.22) 


V    =  C14^^  -  C15;^^  (C.23) 

and 

-  PeClbC^j;  0^  +  ij;  i  )  +  COe^  (C.24) 

Rearranging  the  right-hand  sides  of  Eqs.  (C.20)  and  (C.24) 
yield 

[C7  +  RaiCM^,   -   C18^  )]q.     . 

=  (C8  -  ReClS^;^)^.  .    +  (C9  +  ReCl6i  )n.  ,  , 

+  (CIO  +  ReC^bi,^)^^^^^^    +    (C11  -  ReCIbi  )fii  ^_-i 

+  C05^^  (C.25) 

and 

Cle.^  .  =  (C8  -  PeCl6;^^)o.,,^j  .  (C9  .   PeCl6^^ )  e^.,  ^  . 

+  (C10  +  PeCl6;i;  )0      +  (C11  -  PeClbij;  )0.  ,  , 

+  COi^^  (C.26) 

Equations  (C.21)  to  (C.23)  and  (C.25)  and  (C.2b)  are  the 
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central-difference  versions  of  the  original  governing 
partial -differential-equations. 

The  numerical  solution  of  Eqs.  (C.25)  and  (C.26)  may- 
present  a  problem  because  it  does  not  always  converge.   In 
order  to  rectify  this  problem,  one  way  to  do  it  is  to  make 
the  diagonal  dominant  in  the  matrix  expressions  for  these 
equations.   Notice  that  the  up-wind  method  cannot  be 
directly  used  here  because  the  convection  terms  in  those 
transport  equations  have  been  simplified  to  the  extent  that 
their  original  formats  have  lost.   Use  is  thus  made  of 
Takemitsu's  method  which  is  slightly  modified  as  presented 
below. 

There  are  four  terms  on  the  right-hand  side  of  Eq. 
(C.25)  that  can  be  regrouped  to  relate  convection  as 


CV  =  Re[  -  C^6^^{n..     .    -   ^.    .     ■) 


+  Cl6^,(n.  .  .  -  fi.  .  .)]  (C.27) 


Introducing 


U  =  C16.)^   ,  V  =  -  C16^^  (C.28) 


permits  rewriting  Eq.  (C.27)  more  compactly  as 

CV  ,  Ret-  U(«.^,_.  -  a..^_j)  -  V(a._.^^  -  a._._,)]       (c.29) 
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Notice  that  C16  is  always  positive,  while  ^      and  ^^    may 
change  sign,  depending  on  the  local  stream  function.   The 
terms  in  Eq.  (C.29)  have  been  written  in  central  difference 
forms  and  they  must  be  changed  to  either  forward  or 
backward  difference  to  assure  a  convergent  solution. 
Experience  shows  its  choice  hinges  on  the  sign  of  U  and  V; 
a  forward  difference  must  be  used  if  either  of  them  is 
negative  and  vice  versa.   To  meet  these  requirements,  Eq. 
(C.29)  is  written  as 


i,j  +  1     i,j     i,j     i,j-V-' 


(C.30) 


which  reduces  to  the  following  four  situations 
(i)    U  and  V  are  both  positive: 


CV  =  2Re[-  \]{n.     .    -    n.    ^     .)  -  V(n.  .  -  n.  .  J]   (C.31) 
iij     1-1, J       i,J     i,J-1 


(ii)    U  is  positive  and  V  is  negative 


CV  =  2Re[-  U(n.  .  -  n.  ,  .)  -  V(f^.  .  .  -  n.  .)]   (C.32) 


(iii)   U  is  negative  and  V  is  positive 


CV    =   2Re[-   \]{a.     .     .    -    n.     .)    -   V(n.     .    -    ^.     .    ,)]       (C.33) 
1+1,0  1,0  1,0  1,0-1 
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(iv)        U   and   V   are    both  negative: 

CV    =   2Re[-   U(fi.    .    .    -    ^.     .)    -   V(j^.     .    .    -    n.     .)]       (C.34) 
1+1,  J  ijj  -i-»J"'"'  ■'■to 

The  general  form  of  CV,  Eq.  (C.30),  is  now  substituted 
in  Eq.  (C.25)  and  with  some  regrouping  of  terms  changes  the 
latter    to 

{C7    +   Re[C17?,   -   C184i      +    2(|u|    +    |v|)]}q.     . 
^  5  n  1     I         I      I  1  >  0 


=    [C8   -   Re(U   -      U    )]a.    ,     .    +    [C9   +   Re(U    +      U    ) ] n .    .     . 
II         1+1, J  II         -'-~i,j 


+    [CIO    -    Re(V    -      V    )]n.     .    .    +    [C11    +    Re(V    + 
I       I         1 ,  J  +  I 

+    CO^, 


V  hjfi 


5ti 


i,d-1 
(C.35) 


Following  the  similar  approach,  a  general  form  of  the 
energy  equation  is  derived  as 


[CI    +    2Pe(    U      +      V    ) ]0.     . 
II         II         1 ,  J 


=    [C8   -   Pe(U    -      u|)]0.    ,    ,    +    [C9   +   Pe(U    +      U  hje 

|l  1+1,  J  II  -LI,  J 

+    [C10   -   Pe(V   -  I  v|)]0.     .    ,    +    [C11    +   Pe    (V    +    |v|)]0.     .    . 
II         i,j+'  i,j~i 


+    CO0 


5n 


(C.36) 


These  equations  will  lead  to  a  convergent  solution. 

In  closing  this  appendix,  a  short  note  is  included  to 
show  that  the  scneme  developed  above  is  rooted  in  the  up- 
wind scheme  mentioned  earlier.   According  to  this  scheme. 
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it  is  necessary  to  consider  the  flow  retains  its  past 
identity  when  it  is  migrated  to  a  new  location.   The 
backward  and  forward  represei^tations  of  ^   discussed  above 
indeed  come  from  this  origin.   In  fact,  it  can  be  shown 
that  the  sign  of  U  and  V  are  related  to  that  of  u  and  v, 
which  can  be  verified  by  using  those  equations  in  domain  B, 


APPENDIX  D 
DERIVATION  OF  PRESSURE  GRADIENT  EQUATIONS 


The  expressions  for  pressure  gradients  in  z  and  r 
directions  at  the  wall  will  be  derived  in  this  appendix. 
In  the  derivation,  the  Navier-Stokes  equations  will  be  used 
together  with  the  continuity  equation  to  derive  the 
pressure  gradient  expressions.   Finally,  the  vorticity 
definition  is  introduced  to  reduce  the  differential  order 
in  the  final  expressions. 

For  axisyraraetrical  flow  of  an  incompressible  medium  in 
a  cylindrical  pipe  or  annulus,  the  nondimensional 
continuity  equation  and  momentum  equations  are  given  as 
follows: 

8z    3r    r 

3u      au      9p     1  , 9^u    9^u    1  9u.  ,^    ^. 

and 

9v      9v      3d     1  ,9^v    9^v    1  9v     v^        .^  ^. 
dz  3r      3r    Re   g^2    g^2    r  3r    ^2 


where 


■X-      ^o 

P  =  P  /(P^m  ) 
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Equation  (D.2)  is  used  to  write  the  pressure  gradient 
In  z  direction  at  the  wall  as 

|L.   1  (!%,!%.  1|^)  (D.4) 

The  continuity  equation  is  differentiated  with  respect 
to  z 


8^U       8  V    1  3v 


^  2      3z9r    r  3z 

dZ 


and  substituted  into  (D.4)  to  give 


(D.5) 


3p 1  ( 9^u    a^v   J_  in  _  1  iv_^  n^  K ^ 

3z  ~  Re  \  2  "  az8r  "^  r  3r    r  ^z^  Ku.o) 

dr 


Recall  that  the  vorticity  is  defined  as 


=  =lf-17  (D-7) 


Differentiating  this  equation  with  respect  to  r  gives 


ar  -  3z3r  "  ,  2  ^'^•^^ 

3r 


Introducing  these  two  equations  into  (D.6)  gives  the  final 
expression  for  the  pressure  gradient  in  z  direction  at  the 
wall 


121 


3z     Re  ^r   ar^  ^^•y'' 


The  pressure  gradient  in  r  direction  at  the  wall  may 
be  derived  by  following  a  similar  procedure.   Equation 
(D.3)  will  be  used  instead,  and  this  equation  is  simplified 
by  using  the  condition  that  u  =  v  =  0  at  the  wall.   It 
follows  that 


3r   Re  \  2     2  ^  r  ar-"  ^^*  '^^ 

8z    3r 


Equation  (D.I)  is  differentiated  with  respect  to  r 


a^u   a^v^ii^.v  ^^^^^^ 


ar9z    g^2    r  ar    ^2 


and  Eq.  (D.7)  is  differentiated  with  respect  to  z  to  give 


2      2 
in  _  9  V    a  u  ,  . 

az  -  g^2  "  araz  ^^- '^^ 


and  finally  introducing  these  two  equations  into  (D.10) 
gives  the  pressure  gradient  in  r  direction  at  the  wall  as 


i£.  _  _L  iiL  fn  1^ 

a  r    Re  a  z  ^  u .  I :) 
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